#
# Proj4.pd
#
# PDL::PP Definition file for the PDL::Transform::Proj4 module
#
# Judd Taylor, USF IMaRS
# 4 Apr 2006
#

use strict;

# It won't build without this:
# 8 Aug 2006 - JT - Well, I guess it can!
#use lib '../../../blib/lib';
#use lib '../../../blib/arch';

use PDL;
#use PDL::Transform; # This is not needed, and just causes trouble!
use PDL::GIS::Proj;

use vars qw( $VERSION );
$VERSION = "1.32";

pp_add_isa( 'PDL::Transform' );

# This array holds the list of functions we want to export (everything here is explicit!)
#
my @export_funcs = ();

pp_addbegin( <<'ENDBEGIN' );
use PDL;
use PDL::NiceSlice;
use PDL::Transform;
use PDL::GIS::Proj;
ENDBEGIN

#
# Put in the general projection:
#
pp_addpm( { At => 'Top' }, <<'ENDPM' );

#
# PDL::Transform::Proj4
#
# Judd Taylor, USF IMaRS
# 4 Apr 2006
#

=head1 NAME

PDL::Transform::Proj4 - PDL::Transform interface to the Proj4 projection library

=head1 SYNOPSIS

 # Using the generalized proj interface:
 # Make an orthographic map of Earth
 use PDL::Transform::Cartography;
 use PDL::Transform::Proj4;
 $a = earth_coast();
 $a = graticule(10,2)->glue(1,$a);
 $t = t_proj( proj_params => "+proj=ortho +ellps=WGS84 +lon_0=-90 +lat_0=40" );
 $w = pgwin(xs);
 $w->lines($t->apply($a)->clean_lines());
 
 # Using the aliased functions:
 # Make an orthographic map of Earth
 use PDL::Transform::Cartography;
 use PDL::Transform::Proj4;
 $a = earth_coast();
 $a = graticule(10,2)->glue(1,$a);
 $t = t_proj_ortho( ellps => 'WGS84', lon_0 => -90, lat_0 => 40 )
 $w = pgwin(xs);
 $w->lines($t->apply($a)->clean_lines());

=head1 DESCRIPTION

Works like PDL::Transform::Cartography, but using the proj library in the background.

Please see the proj library docs at L<http://www.remotesensing.org/proj> for more information
on proj, and how to use the library.

=head1 GENERALIZED INTERFACE

The main object here is the PDL::Transform::Proj4 object, aliased to the t_proj() function.

This object accepts all of the standard options described below, but mainly is there to be called
with just the B<proj_params> option defined.

When options are used, they must be used with a '+' before them when placed in the proj_params string,
but that is not required otherwise. See the SYNOPSIS above.

=head2 ALIASED INTERFACE

Other than t_proj(), all of the other transforms below have been autogenerated, and may not work
properly. The main problem is determining the parameters a projection requires from the proj
library itself. 

Due to the difficulties in doing this, there may be times when the proj docs specify a parameter 
for a projection that won't work using the anon-hash type specification. In that case, just throw
that parameter in the proj_params string, and everything should work fine.

=head1 PARAMETERS AVAILABLE IN ALL PROJECTIONS

=head2 General Parameters

=head3 proj_params

This is a string containing the proj "plus style" parameters. This would be similar to what you 
would put on the command line for the 'proj' tool. Like "+proj=ortho +ellps=WGS84 +lon_0=-90 +lat_0=40".

This parameter overrides the others below when it contains parameters that are also specified 
explicitly.

=head3 proj

The proj projection code to use (like ortho...)

=head3 x_0

Cartesian X offset for the output of the transformation

=head3 y_0

Cartesian Y offset for the output of the transformation

=head3 lat_0

Central latitude for the projection.
NOTE: This may mean other things depending on the projection selected, read the proj docs!

=head3 lon_0

Central longitude for the projection.
NOTE: This may mean other things depending on the projection selected, read the proj docs!

=head3 units

Cartesian units used for the output of the projection.
NOTE: Like most of the options here, this is likely useless in the current implementation
of this library.

=head3 init

Specify a file:unit for proj to use for its runtime defaults. See the proj docs.

=head3 no_defs

Don't load any defaults. See the proj docs.

=head3 over

Normally, the transformation limits the output to between -180 and 180 degrees (or the
cartesian equivalent), but with this option that behavior is turned off.

=head3 geoc

Input values are geocentric coordinates.

=head2 Earth Figure Parameters

=head3 ellps

Ellipsoid datum to use. Ex: WGS72, WGS74.
See the proj docs and command line tool for list of possibilities ('proj -le').

=head3 R

Radius of the Earth.

=head3 R_A

Radius of a sphere with equivalent surface area of specified ellipse. 

=head3 R_V

Radius of a sphere with equivalent volume of specified ellipse. 

=head3 R_a

Arithmetic mean of the major and minor axis, Ra = (a + b)/2. 

=head3 R_g

Geometric mean of the major and minor axis, Rg = (ab)1/2. 

=head3 R_h

Harmonic mean of the major and minor axis, Rh = 2ab/(a + b). 

=head3 R_lat_a=phi

Arithmetic mean of the principle radii at latitude phi.

=head3 R_lat_g=phi

Geometric mean of the principle radii at latitude phi.

=head3 b

Semiminor axis or polar radius 

=head3 f

Flattening

=head3 rf

Reciprocal flattening, +rf=1/f

=head3 e

Eccentricity +e=e 

=head3 es

Eccentricity squared +es=e2

=cut

sub new 
{
    my $proto = shift;
    my $sub = "PDL::Transform::Proj4::new()";
    #print STDERR "$sub: ARGS: [" . join(", ", @_ ) . "]\n";
    my $class = ref($proto) || $proto;
    my $self  = $class->SUPER::new( @_ );
    
    bless ($self, $class);
    
    my $o = $_[0];
    unless( (ref $o) ) 
        { $o = {@_}; }
    
    #use Data::Dumper;
    #my $dd2 = Data::Dumper->new( [$o], ["$sub: o"] );
    #$dd2->Indent(1);
    #print STDERR $dd2->Dump();
        
    $self->{name} = "Proj4";

    # Grab our options:
    
    # Used in the general sense:
    $self->{params}->{proj_params} = PDL::Transform::_opt( $o, ['proj_params','params'] );

    # Projection options available to all projections:
    $self->{general_params} = [ qw( proj x_0 y_0 lat_0 lon_0 units init ) ];
    foreach my $param ( @{ $self->{general_params} } )
        { $self->{params}->{$param} = PDL::Transform::_opt( $o, [ $param ] ); }
    
    # Options that have no value (like "+over"):
    $self->{bool_params} = [ qw( no_defs over geoc ) ];
    foreach my $param ( @{ $self->{bool_params} } )
        { $self->{params}->{$param} = ( PDL::Transform::_opt( $o, [ $param ] ) ) ? 'ON' : undef; }
    
    # Options for the Earth figure: (ellipsoid, etc):
    $self->{earth_params} = [ qw( ellps R R_A R_V R_a R_g R_h R_lat_a R_lat_g b f rf e es ) ];
    foreach my $param ( @{ $self->{earth_params} } )
        { $self->{params}->{$param} = PDL::Transform::_opt( $o, [ $param ] ); }
    
    # First process the old params that may already be in the string:
    # These override the specific params set above:
    if( defined( $self->{params}->{proj_params} ) )
    {
        $self->{orig_proj_params} = $self->{params}->{proj_params};
        
        my @params = split( /\s+/, $self->{orig_proj_params} );
        foreach my $param ( @params )
        {
            if( $param =~ /^\+(\S+)=(\S+)/ )
            {
                my ($name, $val) = ($1, $2);
                $self->{params}->{$name} = $val;
                #print STDERR "$sub: $name => $val\n";
            }
            elsif( $param =~ /^\+(\S+)/ )
            {   # Boolean option
                $self->{params}->{$1} = 'ON';
            }                
        }
    }
    
    # Update the proj_string to current options:
    #
    $self->update_proj_string();
    
    #my $dd = Data::Dumper->new( [$self->{params}], ["$sub: params"] );
    #$dd->Indent(1);
    #print STDERR $dd->Dump();
    
    ##############################
    # The meat -- just copy and paste from Transform.pm :)
    #    (and do some proj stuff here as well)
    
    # Forward transformation:
    $self->{func} = sub 
    {
        my $in = shift;
        my $opt = shift;
        my $sub = "PDL::Transform::Proj4->{func}()";
        
        my $out = $in->new_or_inplace();
        
        # Always set the badflag to 1 here, to handle possible bad projection values:
        $out->badflag(1);
        
        PDL::GIS::Proj::fwd_trans_inplace( $out->((0)), $out->((1)), $opt->{proj_params}, 1 );
        return $out;
    };
  
    # Inverse transformation:
    $self->{inv} = sub 
    {
        my $in = shift;
        my $opt = shift;
        my $sub = "PDL::Transform::Proj4->{inv}()";
        
        my $out = $in->new_or_inplace();
        
        # Always set the badflag to 1 here, to handle possible bad projection values:
        $out->badflag(1);
        
        PDL::GIS::Proj::inv_trans_inplace( $out->((0)), $out->((1)), $opt->{proj_params}, 1 );
        return $out;
    };
  
    return $self;
} # End of new()...

sub update_proj_string
{
    my $self = shift;
    my $sub = "PDL::Transform::Proj4::update_proj_string()";
    
    # (Re)Generate the proj_params string from the options passed:
    #
    delete( $self->{params}->{proj_params} );
    my $proj_string = "";
    
    foreach my $param ( keys %{ $self->{params} } )
    {
        next unless defined( $self->{params}->{$param} );
        
        $proj_string .= ( $self->{params}->{$param} eq 'ON' ) 
                        ? "+$param " : " +$param=" . $self->{params}->{$param} . " ";
        #print STDERR "$sub: Adding \'$proj_string\'...\n";
    }
    
    #print STDERR "$sub: Final proj_params: \'$proj_string\'\n";
    
    $self->{params}->{proj_params} = $proj_string;
} # End of update_proj_string()...

sub proj_params
{
    my $self = shift;
    $self->update_proj_string();
    return $self->{params}->{proj_params};
} # End of proj_params()...

sub t_proj 
{ 
    PDL::Transform::Proj4->new( @_ );
} # End of t_proj()...

1;

ENDPM

#
# Add the docs for t_proj:
#
pp_addpm( { At => 'Middle' }, <<'ENDPM' );

=head1 FUNCTIONS

=head2 t_proj

This is the main entry point for the generalized interface. See above on its usage.

=cut
ENDPM
push( @export_funcs, 't_proj' );

#
# Add in the auto-getnerated projection classes:
#
my $projections = PDL::GIS::Proj::load_projection_information();

foreach my $name ( sort keys %$projections )
{
    #print STDERR "Generating code for projection $name...\n";
    
    my $projection = $projections->{$name};

    # Start out with  a blank template:
    my $template = <<'ENDTEMPLATE';

#
# Autogenerated code for the Proj4 projection code: 
#    INSERT_NAME_HERE
#
package PDL::Transform::Proj4::INSERT_NAME_HERE;
use PDL::Transform::Proj4;
@ISA = ( 'PDL::Transform::Proj4' );

sub new
{
    my $proto = shift;
    my $class = ref($proto) || $proto;
    my $sub = "PDL::Transform::Proj4::INSERT_NAME_HERE::new()";
    #print STDERR "$sub: ARGS: [" . join(", ", @_ ) . "]\n";
    my $self  = $class->SUPER::new( @_ );
    bless ($self, $class);
    
    my $o = $_[0];
    unless( (ref $o) ) 
        { $o = {@_}; }
    
    #use Data::Dumper;
    #my $dd2 = Data::Dumper->new( [$o], ["$sub: o"] );
    #$dd2->Indent(1);
    #print STDERR $dd2->Dump();
        
    $self->{name} = "INSERT_FULL_NAME_HERE";
    $self->{proj_code} = "INSERT_NAME_HERE";

    # Make sure proj is set in the options:
    $self->{params}->{proj} = $self->{proj_code};
    
    # Grab our projection specific options:
    #
    $self->{projection_params} = [ qw( INSERT_PARAM_LIST_HERE ) ];
    foreach my $param ( @{ $self->{projection_params} } )
        { $self->{params}->{$param} = PDL::Transform::_opt( $o, [ $param ] ); }
    
    $self->update_proj_string();
    
    #my $dd = Data::Dumper->new( [$self->{params}], ["$sub: params"] );
    #$dd->Indent(1);
    #print STDERR $dd->Dump();
    
    #print STDERR "$sub: Final proj_params: \'" . $self->{params}->{proj_params} . "\'\n";
        
    return $self;
} # End of PDL::Transform::INSERT_NAME_HERE::new()...

1;

ENDTEMPLATE

    # Fill in the projection name:
    $template =~ s/INSERT_NAME_HERE/$name/sg;
    
    # Fill in the full name of the projection:
    $template =~ s/INSERT_FULL_NAME_HERE/$projection->{NAME}/sg;
    
    # Fill in the parameter list:
    my $param_list = join(' ', @{ $projection->{PARAMS}->{PROJ} } );
    $template =~ s/INSERT_PARAM_LIST_HERE/$param_list/sg;

    # Add the code to the module:
    pp_addpm( {At => 'Bot'}, $template );
    
    # Generate the alias sub:
    my $alias_name = "t_proj_$name";
    push( @export_funcs, $alias_name );
    
    
    my $doc_param_list = "";
    if( scalar( @{ $projection->{PARAMS}->{PROJ} } ) )
    {
        $doc_param_list .= "\n=head3 Projection Parameters\n\n";
        foreach my $param ( sort @{ $projection->{PARAMS}->{PROJ} } )
            { $doc_param_list .= "=head4 $param\n\n"; }
    }
    
    my $alias_template = <<'ENDTEMPLATE';
=head2 INSERT_ALIAS_NAME_HERE

Autogenerated transformation function for Proj4 projection code INSERT_NAME_HERE. 

The full name for this projection is INSERT_FULL_NAME_HERE.
INSERT_PARAM_LIST_HERE
=cut
sub INSERT_ALIAS_NAME_HERE 
    { PDL::Transform::Proj4::INSERT_NAME_HERE->new( @_ ); }
ENDTEMPLATE
    
    $alias_template =~ s/INSERT_ALIAS_NAME_HERE/$alias_name/sg;
    $alias_template =~ s/INSERT_NAME_HERE/$name/sg;
    $alias_template =~ s/INSERT_FULL_NAME_HERE/$projection->{NAME}/sg;
    $alias_template =~ s/INSERT_PARAM_LIST_HERE/$doc_param_list/sg;

    pp_addpm( {At => 'Middle'},  $alias_template );

} # for each projection...

pp_add_exported('', join(' ', @export_funcs ) );


# Empty pp_def(), so it will actually generate the code below in the output file!
#
pp_def( '_proj4_dummy',
        Pars => 'i(); [o] o();',
        Doc => undef,
        Code => ';' );

# Add the end docs:
#
pp_addpm( {At => 'Bot'}, <<'ENDDOC');

=head1 AUTHOR & MAINTAINER

Judd Taylor, Orbital Systems, Ltd.
judd dot t at orbitalsystems dot com

=cut

ENDDOC
        
pp_done();