#
# Proj.pd - PP def file for the PROJ->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 warnings;
our $VERSION = "1.32";
pp_addbegin (<<'EOAP');
use Alien::proj;
if ($^O =~ /MSWin32/ and $Alien::proj::VERSION le '1.25') {
$ENV{PATH} = join ';', (Alien::proj->bin_dirs, $ENV{PATH});
}
EOAP
pp_addpm({At=>'Top'},<<'EODOC');
use strict;
use warnings;
=head1 NAME
PDL::GIS::Proj - PDL interface to the PROJ projection library.
=head1 DESCRIPTION
For more information on the PROJ library, see: L<http://www.proj.org/>
EODOC
pp_addpm({At=>'Bot'},<<'EODOC');
=head1 AUTHOR
Judd Taylor, Orbital Systems, Ltd.
judd dot t at orbitalsystems dot com
=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.
=cut
EODOC
#
# Header files:
#
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
pp_addpm(<<'EOPM');
=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('', ' get_proj_info ');
my $code_head = <<'EOF';
PJ_COORD in, out;
PJ *proj = proj_create( NULL, $COMP(params) ); /* Init the projection */
if( proj == NULL )
$CROAK("%s: Projection initialization failed: %s\n", func, proj_errno_string(proj_errno(proj)));
/* Loop over the values converting as we go */
broadcastloop %{
EOF
my $code_foot = <<'EOF';
%}
proj_destroy(proj);
EOF
sub wrap_code {
my ($name, $in, $out, $is_fwd) = @_;
my $body = <<EOF;
PDL_IF_BAD(
if ( @{[ join '||', map "\$ISBAD($in(n=>$_))", 0,1 ]} )
{
@{[ join "\n", map "\$SETBAD($out(n=>$_));", 0,1 ]}
}
else,)
{
in.uv.u = \$$in(n=>0) @{[$is_fwd ? ' * DEG_TO_RAD' : '']};
in.uv.v = \$$in(n=>1) @{[$is_fwd ? ' * DEG_TO_RAD' : '']};
out = proj_trans(proj, PJ_@{[$is_fwd ? 'FWD' : 'INV']}, in);
if (out.uv.u == HUGE_VAL)
{
PDL_IF_BAD(
@{[ join "\n", map "\$SETBAD($out(n=>$_));", 0,1 ]}
,
\$CROAK("%s: Projection conversion failed at (%f, %f): %s\\n",
func, \$$in(n=>0), \$$in(n=>1), proj_errno_string(proj_errno(proj)));
)
}
\$$out(n=>0) = out.uv.u@{[$is_fwd ? '' : ' * RAD_TO_DEG']};
\$$out(n=>1) = out.uv.v@{[$is_fwd ? '' : ' * RAD_TO_DEG']};
}
EOF
my $code = join '', qq{char* func = "$name" PDL_IF_BAD("[BADCODE]",) "()";},
$code_head,
$body,
$code_foot,
;
(Code => $code,
GenericTypes => ['D'],
OtherPars => 'char* params;',
HandleBad => 1,
Inplace => 1,
);
}
#
# Forward transformation:
#
pp_def( 'fwd_transform',
Pars => 'lonlat(n=2); [o] xy(n);',
Doc => <<'EOF',
=for ref
PROJ forward transformation $params is a string of the projection
transformation parameters.
Returns a pdl with x, y values at positions 0, 1. The units are dependent
on PROJ behavior. They will be PDL->null if an error has occurred.
=for bad
Ignores bad elements of $lat and $lon, and sets the corresponding elements
of $x and $y to BAD
EOF
wrap_code('fwd_transform', 'lonlat', 'xy', 1),
);
#
# Inverse Transformation:
#
pp_def( 'inv_transform',
Pars => 'xy(n=2); [o] lonlat(n);',
Doc => <<'EOF',
=for ref
PROJ inverse transformation $params is a string of the projection
transformation parameters.
Returns a pdl with lon, lat values at positions 0, 1. The units are
dependent on PROJ behavior. They will be PDL->null if an error has
occurred.
=for bad
Ignores bad elements of $lat and $lon, and sets the corresponding elements
of $x and $y to BAD
EOF
wrap_code('inv_transform', 'xy', 'lonlat', 0),
);
#
# Utility functions for getting projection description information (in a general case).
#
pp_addxs('', <<'ENDXS' );
HV*
load_projection_descriptions()
CODE:
const PJ_OPERATIONS *lp;
SV* scalar_val;
RETVAL = newHV();
for (lp = proj_list_operations() ; lp->id ; ++lp)
{
scalar_val = newSVpv( *lp->descr, 0 );
hv_store( RETVAL, lp->id, strlen( lp->id ), scalar_val, 0 );
}
OUTPUT:
RETVAL
void
proj_version(...)
PPCODE:
EXTEND(sp, 3);
mPUSHu(PROJ_VERSION_MAJOR);
mPUSHu(PROJ_VERSION_MINOR);
mPUSHu(PROJ_VERSION_PATCH);
ENDXS
pp_add_exported('', ' load_projection_descriptions proj_version ');
pp_addpm( <<'ENDPM' );
=head2 proj_version
Returns a 3-element list with PROJ major, minor, patch version-numbers.
=cut
my %SKIP = map +($_=>1), qw(
and or Special for Madagascar
fixed Earth For CH1903
);
sub load_projection_information
{
my $descriptions = PDL::GIS::Proj::load_projection_descriptions();
my $info = {};
foreach my $projection ( sort keys %$descriptions )
{
my $description = $descriptions->{$projection};
my $hash = {CODE => $projection};
my @lines = split( /\n/, $description );
chomp @lines;
# 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:
$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:
$hash->{PARAMS}{PROJ} = [
grep !$SKIP{$_}, map {s/=//; s/[,\[\]]//sg; $_}
grep length, map split(/\s+/), @lines
];
# 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();