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

use strict;
use warnings;
use File::Spec;
use IPC::Run qw(run);

our $VERSION = "2.098";
pp_setversion($VERSION);

pp_add_isa( 'PDL::Transform' );

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

pp_addhdr(<<'EOHDR');
#include "proj.h"
#include <string.h>
/* from proj_api.h */
#define RAD_TO_DEG      57.295779513082321
#define DEG_TO_RAD      .017453292519943296
EOHDR

sub wrap_code {
  my ($name, $in, $out, $is_fwd) = @_;
  pp_def($name,
    Code => <<EOF,
PJ_CONTEXT *pj_context = proj_context_create();
if (!pj_context)
  \$CROAK("Projection context creation failed: %s\\n", proj_errno_string(proj_errno(NULL)));
PJ *proj = proj_create( pj_context, \$COMP(params) );
if (!proj)
  \$CROAK("Projection initialization failed: %s\\n", proj_errno_string(proj_errno(NULL)));
char rad_in = @{[$is_fwd ? 'proj_angular_input(proj, PJ_FWD)' : '0']},
  deg_out = @{[$is_fwd ? '0' : 'proj_angular_output(proj, PJ_INV)']};
loop (n) %{
  PJ_COORD in, out;
  PDL_IF_BAD(
  if ( @{[ join '||', map "\$ISBAD($in(ncoord=>$_))", 0,1 ]} ) {
    loop (ncoord) %{ \$SETBAD($out()); %}
    continue;
  },)
  loop (ncoord) %{ in.v[ncoord] = rad_in ? DEG_TO_RAD*\$$in() : \$$in(); %}
  out = proj_trans(proj, PJ_@{[$is_fwd ? 'FWD' : 'INV']}, in);
  if (out.v[0] == HUGE_VAL) {
    PDL_IF_BAD(
    loop (ncoord) %{ \$SETBAD($out()); %}
    continue;
    ,
    \$CROAK("Projection conversion failed at (%f, %f): %s\\n",
      \$$in(ncoord=>0), \$$in(ncoord=>1), proj_errno_string(proj_errno(proj)));
    )
  }
  loop (ncoord) %{ \$$out() = deg_out ? RAD_TO_DEG*out.v[ncoord] : out.v[ncoord]; %}
%}
proj_destroy(proj);
proj_context_destroy(pj_context);
EOF
    Pars => "$in(ncoord=2,n); [o] $out(ncoord,n);",
    GenericTypes => ['D'],
    OtherPars => 'char* params;',
    HandleBad => 1,
    Inplace => 1,
    Doc => undef,
  );
}

wrap_code('fwd_transform', 'lonlat', 'xy', 1);
wrap_code('inv_transform', 'xy', 'lonlat', 0);

# Utility functions for getting projection description information (in a general case).
pp_addxs('', <<'ENDXS' );
void
proj_version(...)
  PPCODE:
    EXTEND(sp, 3);
    mPUSHu(PROJ_VERSION_MAJOR);
    mPUSHu(PROJ_VERSION_MINOR);
    mPUSHu(PROJ_VERSION_PATCH);

# returns input_units, output_units
void
units(proj_str)
  char *proj_str
PPCODE:
  PJ *proj = proj_create( NULL, proj_str ); /* Init the projection */
  if (!proj)
    croak("Failed to create PJ from '%s': %s", proj_str, proj_errno_string(proj_errno(proj)));
  EXTEND(sp, 2);
  char *input_u = proj_angular_input(proj, PJ_FWD) ||
    proj_degree_input(proj, PJ_FWD) ? "degrees" : "metres";
  char *output_u = proj_angular_output(proj, PJ_FWD) ||
    proj_degree_output(proj, PJ_FWD) ? "degrees" : "metres";
  mPUSHs(newSVpv(input_u, 0));
  mPUSHs(newSVpv(output_u, 0));
  proj_destroy(proj);
ENDXS

my %SKIP = map +($_=>1), qw(
  and or Special for Madagascar
  fixed Earth For CH1903
);
sub load_projection_information {
    my @cmd = ("proj", "-lP");
    my ($text, $stderr, $exit_code, $in);
    run \@cmd, \$in, \$text, \$stderr;
    $exit_code = $?;
    warn $stderr if $stderr;
    die "proj -lP error $exit_code. See above for error text." if $exit_code;
    my @chunks = $text =~ /(.+?)(?=(?:^\S|\z))/gms;
    chomp for @chunks;
    my %descriptions = map {
      my ($id, $rest) = split /\s*:\s*/, $_, 2;
    } @chunks;
    my %info;
    foreach my $projection ( sort keys %descriptions ) {
        my $description = $descriptions{$projection};
        my %hash = (CODE => $projection);
        my @lines = split( /\n/, $description );
        chomp @lines;
        # Can this projection do inverse?
        $hash{INVERSE} = 1;
        # Full name of this projection:
        ($hash{NAME}, my $temp) = splice @lines, 0, 2;
        if ($temp) {
          # The second line is usually a list of projection types this one is:
          $hash{INVERSE} = 0 if $temp =~ s/no inv\.*,*//i;
          $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:
        # Projection Specific Parameters:
        $hash{PARAMS}{PROJ} = [
          grep !$SKIP{$_}, map {s/=//; s/[,\[\]]//sg; $_}
            grep length, map split(/\s+/), @lines
        ];
        $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;
}

pp_addpm(<<'ENDPM');
=head2 proj_version

Returns a 3-element list with PROJ major, minor, patch version-numbers.
Not exported.

=cut

ENDPM
pp_addpm("sub load_projection_information {+" . do {
  use Data::Dumper;
  local $Data::Dumper::Sortkeys = 1;
  local $Data::Dumper::Terse = 1;
  local $Data::Dumper::Indent = 1;
  join '', Dumper load_projection_information();
}."}\n\n");
pp_export_nothing();

pp_addpm( { At => 'Top' }, <<'ENDPM' );
# PDL::Transform::Proj4
#
# Judd Taylor, USF IMaRS
# 4 Apr 2006

use strict;
use warnings;
use PDL::LiteF;
use PDL::Transform;

=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;
 use PDL::Graphics::Simple;
 $e = earth_image('day');
 ($c, $pen) = graticule(10,2)->glue(1,earth_coast())->clean_lines;
 $t = t_proj(proj_params => "+proj=ortho +ellps=WGS84 +lon_0=0 +lat_0=40");
 $w = pgswin();
 $w->plot(with=>'fits', $e->map($t),
   with=>'polylines', clean_lines($c->apply($t), $pen), {j=>1});

 # Using the aliased functions:
 # Make an orthographic map of Earth
 use PDL::Transform::Cartography;
 use PDL::Transform::Proj4;
 use PDL::Graphics::Simple;
 $e = earth_image('day');
 ($c, $pen) = graticule(10,2)->glue(1,earth_coast())->clean_lines;
 $t = t_proj_ortho( ellps => 'WGS84', lon_0=>0, lat_0=>40 );
 $w = pgswin();
 $w->plot(with=>'fits', $e->map($t),
   with=>'polylines', clean_lines($c->apply($t), $pen), {j=>1});

=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.

Please note that unlike PROJ, all angles in these operations are
in degrees. This is correctly (as of PDL 2.094) reflected in the
PDL::Transform subclass objects.

=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)/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

# Projection options available to all projections:
our @GENERAL_PARAMS = qw(proj x_0 y_0 lat_0 lon_0 units init);
# Options for the Earth figure: (ellipsoid, etc):
our @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);
# Options that have no value (like "+over"):
our @BOOL_PARAMS = qw( no_defs over geoc );

# The meat -- just copy and paste from Transform.pm :)
#    (and do some proj stuff here as well)
sub _proj4_fwd {
  my ($in, $opt) = @_;
  my $out = $in->new_or_inplace();
  # Always set the badflag to 1 here, to handle possible bad projection values:
  $out->badflag(1);
  $out->inplace(1);
  fwd_transform( $out, $opt->{proj_params} );
}
sub _proj4_inv {
  my ($in, $opt) = @_;
  my $out = $in->new_or_inplace();
  # Always set the badflag to 1 here, to handle possible bad projection values:
  $out->badflag(1);
  $out->inplace(1);
  inv_transform( $out, $opt->{proj_params} );
}

sub new {
    my $self  = shift->SUPER::new( @_ );
    my $o = $_[0];
    $o = {@_} if !ref $o;
    $self->{name} = "Proj4";
    $self->{params}{proj_params} = PDL::Transform::_opt( $o, ['proj_params','params'] );
    foreach my $param ( @BOOL_PARAMS ) {
        my $val = PDL::Transform::_opt( $o, [ $param ] );
        $self->{params}{$param} = 'ON' if $val;
    }
    foreach my $param (@GENERAL_PARAMS, @EARTH_PARAMS) {
        my $val = PDL::Transform::_opt( $o, [ $param ] );
        $self->{params}{$param} = $val if defined $val;
    }
    # 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;
            } elsif ($param =~ /^\+(\S+)/) {   # Boolean option
                $self->{params}{$1} = 'ON';
            }
        }
    }
    @$self{qw(idim odim)} = (2, 2);
    # Update proj_string to current options:
    $self->update_proj_string if $self->{params}{proj};
    $self->{func} = \&_proj4_fwd; # Forward transformation
    $self->{inv} = \&_proj4_inv; # Inverse transformation
    return $self;
} # End of new()...

sub update_proj_string {
    my $self = shift;
    # (Re)Generate the proj_params string from the options passed:
    my @params = map "+$_".($self->{params}{$_} eq 'ON' ? '' : "=$self->{params}{$_}"),
      sort grep $_ ne 'proj_params' && defined $self->{params}{$_}, keys %{$self->{params}};
    $self->{params}{proj_params} = "@params";
  my ($iunit, $ounit) = PDL::Transform::Proj4::units($self->{params}{proj_params});
  @$self{qw(iunit ounit)} = ([($iunit)x2], [($ounit)x2]);
} # 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()...

sub _make_class {
  my ($name, $full_name, $can_inv, $param_list) = @_;
  eval <<"EOF";
package PDL::Transform::Proj4::${name};
our \@ISA = 'PDL::Transform::Proj4';
our \@PARAM_LIST = $param_list;
sub new {
  my \$self  = shift->SUPER::new( \@_ );
  my \$o = \$_[0];
  \$o = {\@_} if !ref \$o;
  \$self->{name} = \$self->{params}{proj} = q{$name};
  foreach my \$param (\@PARAM_LIST) {
      my \$val = PDL::Transform::_opt( \$o, [ \$param ] );
      \$self->{params}{\$param} = \$val if defined \$val;
  }
  delete \$self->{inv} if !$can_inv;
  \$self->update_proj_string;
  return \$self;
}
EOF
}
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-generated projection classes:
my $projections = eval { load_projection_information() } || {};

my $ALIAS_TEMPLATE = pp_line_numbers __LINE__, <<'ENDTEMPLATE';
=head2 INSERT_ALIAS_NAME_HERE

Proj4 projection code C<INSERT_NAME_HERE>, full name "INSERT_FULL_NAME_HERE".
INSERT_CATEGORIES_LIST_HERE
INSERT_PARAM_LIST_HERE

=cut

sub INSERT_ALIAS_NAME_HERE {'PDL::Transform::Proj4::INSERT_NAME_HERE'->new(@_)}
ENDTEMPLATE

foreach my $name ( sort keys %$projections )
{
    my $projection = $projections->{$name};
    my $param_list = 'qw('.join(' ', @{ $projection->{PARAMS}{PROJ} }).')';
    pp_addpm( {At => 'Bot'}, "_make_class(q{$name}, q{$projection->{NAME}}, $projection->{INVERSE}, q{$param_list});\n" );
    # Generate the alias sub:
    my $alias_name = "t_proj_$name";
    push( @export_funcs, $alias_name );
    $_ = "" for my ($doc_param_list, $cat_list);
    if( scalar( @{ $projection->{PARAMS}{PROJ} } ) ) {
        $doc_param_list .= "\nProjection Parameters\n\n=for options\n\n=over 4\n\n";
        $doc_param_list .= join '', map "=item $_\n\n", sort @{ $projection->{PARAMS}{PROJ} };
        $doc_param_list .= "=back\n\n";
    }
    unshift @{$projection->{CATEGORIES}}, "B<no inverse>"
      if !$projection->{INVERSE};
    $cat_list = "\nCategories: @{$projection->{CATEGORIES}}.\n\n"
      if @{$projection->{CATEGORIES} || []};
    my $alias_template = $ALIAS_TEMPLATE;
    $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;
    $alias_template =~ s/INSERT_CATEGORIES_LIST_HERE/$cat_list/sg;
    pp_addpm( {At => 'Middle'},  $alias_template );
} # for each projection...

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

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

=head1 AUTHOR

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

=cut
ENDDOC

pp_done();