#
# Proj.pd - PP def file for the Proj4->PDL interface.
#
# COPYRIGHT NOTICE:
#
# Copyright 2003 Judd Taylor, USF Institute for Marine Remote Sensing (judd@marine.usf.edu).
#
# Now GPL!
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
#
use strict;
use PDL;

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

pp_addpm(<<'EODOC');
=head1 NAME

PDL::GIS::Proj - PDL interface to the Proj4 projection library.

=head1 DESCRIPTION

PDL interface to the Proj4 projection library.

For more information on the proj library, see: http://www.remotesensing.org/proj/

=head1 AUTHOR

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

=head1 DATE

18 March 2003

=head1 CHANGES

=head2 1.32 (29 March 2006) Judd Taylor

    - Getting ready to merge this into the PDL CVS. 
    
=head2 1.31 (???) Judd Taylor

    - Can't remember what was in that version

=head2 1.30 (16 September 2003) Judd Taylor

    - The get_proj_info() function actually works now.

=head2 1.20 (24 April 2003) Judd Taylor

    - Added get_proj_info().

=head2 1.10 (23 April 2003) Judd Taylor

    - Changed from using the proj_init() type API in projects.h to the
    - proj_init_plus() API in proj_api.h. The old one was not that stable...

=head2 1.00 (18 March 2003) Judd Taylor

    - Initial version

=head1 COPYRIGHT NOTICE

Copyright 2003 Judd Taylor, USF Institute for Marine Remote Sensing (judd@marine.usf.edu).

GPL Now!

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

=head1 SUBROUTINES

=cut
EODOC

#
# Header files:
#
pp_addhdr(<<'EOHDR');

#include "projects.h"
#include "proj_api.h"
#include <string.h>

EOHDR

pp_addpm(<<'EOPM');

=head2 fwd_transform($lon(pdl), $lat(pdl), $params)

Proj4 forward transformation $params is a string of the projection transformation
parameters.

Returns two pdls for x and y values respectively. The units are dependant on Proj4
behavior. They will be PDL->null if an error has occurred.

BadDoc: Ignores bad elements of $lat and $lon, and sets the corresponding elements 
of $x and $y to BAD

=cut
sub fwd_transform
{
    my ($lon, $lat, $params) = @_;
    my $x = null;
    my $y = null;

    #print "Projection transformation parameters: \'$params\'\n";

    _fwd_trans( $lon, $lat, $x, $y, $params );

    return ($x, $y);
} # End of fwd_transform()...

=head2 inv_transform($x(pdl), $y(pdl), $params)

Proj4 inverse transformation $params is a string of the projection transformation
parameters.

Returns two pdls for lat and lon values respectively. The units are dependant on Proj4
behavior. They will be PDL->null if an error has occurred.

BadDoc: Ignores bad elements of $lat and $lon, and sets the corresponding elements 
of $x and $y to BAD

=cut
sub inv_transform
{
    my ($x, $y, $params) = @_;
    my $lon = null;
    my $lat = null;

    #print "Projection transformation parameters: \'$params\'\n";

    _inv_trans( $x, $y, $lon, $lat, $params );
    return ($lon, $lat);
} # End of fwd_transform()...

=head2 get_proj_info($params_string)

Returns a string with information about what parameters proj will
actually use, this includes defaults, and +init=file stuff. It's 
the same as running 'proj -v'. It uses the proj command line, so
it might not work with all shells. I've tested it with bash.

=cut
sub get_proj_info
{
    my $params = shift;
    my @a = split(/\n/, `echo | proj -v $params`);
    pop(@a);
    return join("\n", @a);
} # End of get_proj_info()...

EOPM

pp_add_exported('', ' fwd_transform inv_transform get_proj_info ');

#
# Forward transformation:
#
pp_def( '_fwd_trans',
        Pars => 'lon(n);
                 lat(n);
                 [o] x(n);
                 [o] y(n);',
        GenericTypes => ['D'],
        Inplace => ['lon', 'x', 'lat', 'y'],
        OtherPars => 'char* params;',
        HandleBad => 1,
        Doc => undef,
        #BadDoc => 'Ignores bad elements of $lat and $lon, and sets '
        #          . 'the corresponding elements of $x and $y to BAD',
        Code => <<'EOCODE',
/* vars needed: */
char* func = "_fwd_trans()";
char errstr[255];
projUV in, out;
projPJ proj;
const double d2r = DEG_TO_RAD;

/* Init the projection */
proj = pj_init_plus( $COMP(params) );
if( proj == NULL )
{
    sprintf(errstr, "%s: Projection initialization failed: %s\n", func, pj_strerrno(pj_errno));
    croak(errstr);
}

/* Loop over the values converting as we go */
loop (n)
%{
    in.u = $lon() * d2r;
    in.v = $lat() * d2r;

    out = pj_fwd(in, proj);
    if (out.u == HUGE_VAL)
    {
        sprintf(errstr, "%s: Projection conversion failed at (%f, %f): %s\n",
                func, $lon(), $lat(), pj_strerrno(pj_errno));
        croak(errstr);
    }
    $x() = out.u;
    $y() = out.v;
%}
pj_free(proj);
EOCODE


        BadCode => <<'EOBAD');
/* vars needed: */
char* func = "_fwd_trans()";
char errstr[255];
projUV in, out;
projPJ proj;
const double d2r = DEG_TO_RAD;

/* Init the projection */
proj = pj_init_plus( $COMP(params) );
if( proj == NULL )
{
    sprintf(errstr, "%s: Projection initialization failed: %s\n", func, pj_strerrno(pj_errno));
    croak(errstr);
}

/* Loop over the values converting as we go */
loop (n)
%{
    if ( $ISBAD(lon()) || $ISBAD(lat()) )
    {
        $SETBAD(x());
        $SETBAD(y());
    }
    else
    {
        in.u = $lon() * d2r;
        in.v = $lat() * d2r;

        out = pj_fwd(in, proj);
        if (out.u == HUGE_VAL)
        {
            sprintf(errstr, "%s: Projection conversion failed at (%f, %f): %s\n",
                    func, $lon(), $lat(), pj_strerrno(pj_errno));
            croak(errstr);
        }
        $x() = out.u;
        $y() = out.v;
    }
%}
EOBAD

#
# INPLACE Forward transformation: (Call this one directly)
#


pp_addpm( <<'ENDPM' );
#
# Wrapper sub for _fwd_trans_inplace that sets a default for the quiet variable.
# 
sub fwd_trans_inplace
{
    my $lon = shift;
    my $lat = shift;
    my $params = shift;
    my $quiet = shift || 0;
    
    _fwd_trans_inplace( $lon, $lat, $params, $quiet );
} # End of fwd_trans_inplace()...

ENDPM

pp_add_exported('', 'fwd_trans_inplace');

pp_def( '_fwd_trans_inplace',
        Pars => 'lon();
                 lat();',
        GenericTypes => ['F', 'D'],
        OtherPars => 'char* params; int quiet;',
        HandleBad => 1,
        Doc => undef,
        #BadDoc => 'Ignores bad elements of $lat and $lon, and sets '
        #          . 'the corresponding elements of $x and $y to BAD',
        Code => <<'EOCODE',
/* vars needed: */
char* func = "_fwd_trans_inplace()";
char errstr[255];
projUV in, out;
projPJ proj;
const double d2r = DEG_TO_RAD;

/* Init the projection */
proj = pj_init_plus( $COMP(params) );
if( proj == NULL )
{
    sprintf(errstr, "%s: Projection initialization failed: %s\n", func, pj_strerrno(pj_errno));
    croak(errstr);
}

/* Loop over the values converting as we go */
threadloop
%{
    in.u = $lon() * d2r;
    in.v = $lat() * d2r;

    out = pj_fwd(in, proj);
    if (out.u == HUGE_VAL)
    {
        sprintf(errstr, "%s: Projection conversion failed at (%f, %f): %s\n",
                func, $lon(), $lat(), pj_strerrno(pj_errno));
        croak(errstr);
    }
    $lon() = out.u;
    $lat() = out.v;
%}
pj_free(proj);
EOCODE


        BadCode => <<'EOBAD');
/* vars needed: */
char* func = "_fwd_trans_inplace[BADCODE]()";
char errstr[255];
projUV in, out;
projPJ proj;
const double d2r = DEG_TO_RAD;
int loud = ! $COMP(quiet);

/* Init the projection */
proj = pj_init_plus( $COMP(params) );
if( proj == NULL )
{
    sprintf(errstr, "%s: Projection initialization failed: %s\n", func, pj_strerrno(pj_errno));
    croak(errstr);
}

/* Loop over the values converting as we go */
threadloop
%{
    if ( !($ISBAD(lon()) || $ISBAD(lat())) )
    {
        in.u = $lon() * d2r;
        in.v = $lat() * d2r;

        out = pj_fwd(in, proj);
        if (out.u != HUGE_VAL)
        {
            $lon() = out.u;
            $lat() = out.v;
        }
        else
        {
            $SETBAD( lon() );
            $SETBAD( lat() );
            if( loud )
            {
                sprintf(errstr, "%s: Projection conversion failed at (%f, %f): %s\n",
                        func, $lon(), $lat(), pj_strerrno(pj_errno));
                fprintf( stderr, "%s", errstr );
                fprintf( stderr, "%s: NOTE: Subsequent errors may have occured, but I'm only reporting the first!\n", func );   
            }
        }
    }
%}
pj_free(proj);
EOBAD


#
# Inverse Transformation:
#
pp_def( '_inv_trans',
        Pars => 'x(n);
                 y(n);
                 [o] lon(n);
                 [o] lat(n);',
        GenericTypes => ['D'],
        OtherPars => 'char* params;',
        HandleBad => 1,
        Doc => undef,
        #BadDoc => 'Ignores bad elements of $x and $y, and sets '
        #          . 'the corresponding elements of $lon and $lat to BAD',
        Code => <<'EOCODE',
/* vars needed: */
char errstr[255];
projUV in, out;
projPJ proj;
const double r2d = RAD_TO_DEG;

/* Init the projection */
proj = pj_init_plus( $COMP(params) );
if( proj == NULL )
{
    sprintf(errstr, "Projection initialization failed: %s\n", pj_strerrno(pj_errno));
    croak(errstr);
}

/* Loop over the values converting as we go */
loop (n)
%{
    in.u = $x();
    in.v = $y();

    out = pj_inv(in, proj);
    if (out.u == HUGE_VAL)
    {
        sprintf(errstr, "Projection conversion failed: %s\n", pj_strerrno(pj_errno));
        croak(errstr);
    }

    $lon() = out.u * r2d;
    $lat() = out.v * r2d;
%}

pj_free(proj);
EOCODE


        BadCode => <<'EOBAD' );
/* vars needed: */
char errstr[255];
projUV in, out;
projPJ proj;
const double r2d = RAD_TO_DEG;

/* Init the projection */
proj = pj_init_plus( $COMP(params) );
if( proj == NULL )
{
    sprintf(errstr, "Projection initialization failed: %s\n", pj_strerrno(pj_errno));
    croak(errstr);
}

/* Loop over the values converting as we go */
loop (n)
%{
    if ( $ISBAD(x()) || $ISBAD(y()) )
    {
        $SETBAD(lon());
        $SETBAD(lat());
    }
    else
    {
        in.u = $x();
        in.v = $y();

        out = pj_inv(in, proj);
        if (out.u == HUGE_VAL)
        {
            sprintf(errstr, "Projection conversion failed: %s\n", pj_strerrno(pj_errno));
            croak(errstr);
        }

        $lon() = out.u * r2d;
        $lat() = out.v * r2d;
    }
%}
pj_free(proj);
EOBAD

#
# INPLACE Inverse Transformation: (call it directly)
#

pp_addpm( <<'ENDPM' );
#
# Wrapper sub for _fwd_trans_inplace that sets a default for the quiet variable.
# 
sub inv_trans_inplace
{
    my $lon = shift;
    my $lat = shift;
    my $params = shift;
    my $quiet = shift || 0;
    
    _inv_trans_inplace( $lon, $lat, $params, $quiet );
} # End of fwd_trans_inplace()...

ENDPM

pp_add_exported('', 'inv_trans_inplace');

pp_def( '_inv_trans_inplace',
        Pars => 'x();
                 y();',
        GenericTypes => ['F','D'],
        OtherPars => 'char* params; int quiet;',
        HandleBad => 1,
        Doc => undef,
        #BadDoc => 'Ignores bad elements of $x and $y, and sets '
        #          . 'the corresponding elements of $lon and $lat to BAD',
        Code => <<'EOCODE',
/* vars needed: */
char* func = "_inv_trans_inplace()";
char errstr[255];
projUV in, out;
projPJ proj;
const double r2d = RAD_TO_DEG;

/* Init the projection */
proj = pj_init_plus( $COMP(params) );
if( proj == NULL )
{
    sprintf(errstr, "%s: Projection initialization failed: %s\n", func, pj_strerrno(pj_errno));
    croak(errstr);
}

/* Loop over the values converting as we go */
threadloop
%{
    in.u = $x();
    in.v = $y();

    out = pj_inv(in, proj);
    if (out.u == HUGE_VAL)
    {
        sprintf(errstr, "%s: Projection conversion failed: %s\n", func, pj_strerrno(pj_errno));
        croak(errstr);
    }

    $x() = out.u * r2d;
    $y() = out.v * r2d;
%}
pj_free(proj);
EOCODE


        BadCode => <<'EOBAD' );
/* vars needed: */
char* func = "_inv_trans_inplace[BADCODE]()";
char errstr[255];
projUV in, out;
projPJ proj;
const double r2d = RAD_TO_DEG;
int loud = ! $COMP(quiet);

/* Init the projection */
proj = pj_init_plus( $COMP(params) );
if( proj == NULL )
{
    sprintf(errstr, "%s: Projection initialization failed: %s\n", func, pj_strerrno(pj_errno));
    croak(errstr);
}

/* Loop over the values converting as we go */
threadloop
%{
    if ( ! ($ISBAD(x()) || $ISBAD(y())) )
    {
        in.u = $x();
        in.v = $y();

        out = pj_inv(in, proj);
        if (out.u != HUGE_VAL)
        {
            $x() = out.u * r2d;
            $y() = out.v * r2d;
        }
        else
        {
            $SETBAD( x() );
            $SETBAD( y() );
            if( loud )
            {
                /* Don't croak, just set the output to bad: */
                sprintf(errstr, "%s: Projection conversion failed at (%f, %f): %s\n",
                        func, $x(), $y(), pj_strerrno(pj_errno));
                fprintf( stderr, "%s", errstr );
                fprintf( stderr, "%s: NOTE: Subsequent errors may have occured, but I'm only reporting the first!\n", func );
            }
        }
    }
%}
pj_free(proj);
EOBAD

#
# Utility functions for getting projection description information (in a general case).
#

pp_addxs('', <<'ENDXS' );

HV* 
load_projection_descriptions()
    CODE:
        struct PJ_LIST *lp;
        SV* scalar_val;
        RETVAL = newHV();
        
        for (lp = pj_list ; lp->id ; ++lp) 
        {
            scalar_val  = newSVpv( *lp->descr, 0 );
            hv_store( RETVAL, lp->id, strlen( lp->id ), scalar_val, 0 );
        }
    OUTPUT:
        RETVAL

ENDXS
pp_add_exported('', ' load_projection_descriptions ');


#
# Perl code to finish loading the projetion information by parsing the descriptions:
#

pp_addpm( <<'ENDPM' );

sub load_projection_information
{
    my $descriptions = PDL::GIS::Proj::load_projection_descriptions();
    
    my $info = {};
    
    foreach my $projection ( keys %$descriptions )
    {
        my $description = $descriptions->{$projection};
    
        my $hash = {};
        $hash->{CODE} = $projection;
        
        
        
        my @lines = split( /\n/, $description );
        chomp @lines;
        
        # Full name of this projection:
        $hash->{NAME} = $lines[0];
        
        # The second line is usually a list of projection types this one is:
        my $temp = $lines[1];
        $temp =~ s/no inv\.*,*//;
        $temp =~ s/or//;
        my @temp_types = split(/[,&\s]/, $temp );
        my @types = grep( /.+/, @temp_types );
        $hash->{CATEGORIES} = \@types;
        
        # If there's more than 2 lines, then it usually is a listing of parameters:
        
        # General parameters for all projections:
        $hash->{PARAMS}->{GENERAL} = 
            [ qw( x_0 y_0 lon_0 units init no_defs geoc over ) ];
        
        # Earth Figure Parameters:
        $hash->{PARAMS}->{EARTH} = 
            [ qw( ellps b f rf e es R R_A R_V R_a R_g R_h R_lat_g ) ];
        
        # Projection Specific Parameters:
        my @proj_params = ();
        if( $#lines >= 2 )
        {
            foreach my $i ( 2 .. $#lines )
            {
                my $text = $lines[$i];
                my @temp2 = split( /\s+/, $text );
                my @params = grep( /.+/, @temp2 );
                foreach my $param (@params)
                {
                    $param =~ s/=//;
                    $param =~ s/[\[\]]//sg;
                    next if $param =~ /^and$/;
                    next if $param =~ /^or$/;
                    next if $param =~ /^Special$/;
                    next if $param =~ /^for$/;
                    next if $param =~ /^Madagascar$/;
                    next if $param =~ /^fixed$/;
                    next if $param =~ /^Earth$/;
                    next if $param =~ /^For$/;
                    next if $param =~ /^CH1903$/;
                    push(@proj_params, $param);
                }
            }    
        }
        $hash->{PARAMS}->{PROJ} = \@proj_params;
        
        # Can this projection do inverse?
        $hash->{INVERSE} = ( $description =~ /no inv/ ) ? 0 : 1;
        
        $info->{$projection} = $hash;
    }
    
    # A couple of overrides:
    #
    $info->{ob_tran}->{PARAMS}->{PROJ} = 
        [ 'o_proj', 'o_lat_p', 'o_lon_p', 'o_alpha', 'o_lon_c', 
          'o_lat_c', 'o_lon_1', 'o_lat_1', 'o_lon_2', 'o_lat_2' ];
          
    $info->{nzmg}->{CATEGORIES} = [ 'fixed Earth' ];

    return $info;
} # End of load_projection_information()...


ENDPM
pp_add_exported('', ' load_projection_information ');

pp_done();