view release on metacpan or  search on metacpan

lib/PDL/Transform/Proj4.pd  view on Meta::CPAN

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 ($text, $stderr, $exit_code) = Alien::proj->run_utility ("proj", "-lP");
    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;
}