# PDL::PP Definition file for the PDL::Transform::Proj4 module
#
# Judd Taylor, USF IMaRS
# 4 Apr 2006
use strict;
use File::Spec;
use lib File::Spec->catdir((File::Spec->updir) x 3, 'inc');
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-generated projection classes:
require Alien::Proj4;
require PDL::Config;
my @inc = Alien::Proj4->default_inc;
@inc = @{$PDL::Config{PROJ_INC}}
if $PDL::Config{PROJ_INC} and @{$PDL::Config{PROJ_INC}};
my $supplied = File::Spec->catdir((File::Spec->updir) x 2, qw(GIS Proj include));
push @inc, File::Spec->rel2abs($supplied); # because Inline builds elsewhere
Alien::Proj4->import($PDL::Config{PROJ_LIBS}, \@inc);
my $projections = Alien::Proj4->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 .= "\nProjection Parameters\n\n=for options\n\n=over 4\n\n";
foreach my $param ( sort @{ $projection->{PARAMS}->{PROJ} } )
{ $doc_param_list .= "=item $param\n\n"; }
$doc_param_list .= "=back\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();