use Config;
pp_add_isa('PDL::Graphics::PLplot');
pp_addpm({At => Top}, <<'EOD');
$VERSION = '0.10'; # these versions must match!
BEGIN { $VERSION = '0.10'; }; # put this is so it does not use version inherited from PDL::Graphics::PLplot!
=head1 NAME
PDL::Graphics::PLplot::Map - Interface to the GMT coastline database for plotting maps
Perl/PDL interface to GMT's pscoast function to get binary coastline/river/boundary
info into a (big) PDL of latitude/longitude line segments.
NOTE: This module *requires* bad value support! (Use a recent version of PDL.
Set WITH_BADVAL => 1 in perldl.conf).
=head1 SYNOPSIS
use PDL;
use PDL::Graphics::PLplot::Map;
#
## plot just the coastlines in a linear projection
#
my $pl = PDL::Graphics::PLplot::Map->new (DEV => "png", FILE => "map.png");
$pl->worldmap (MAPBOX => [-180, 180, -90, 90]);
#
## Plot coastlines, lat/lon grids and labels in an Azimuthal Equidistant projection
#
$pl = PDL::Graphics::PLplot::Map->new (DEV => "png", FILE => "map1.png");
$pl->worldmap (PROJECTION => 'AZEQDIST', # Azimuthal Equidistant projection
CENTER => [-170, 70], # map centered at 170 deg west lon, 70 deg north lat
RADIUS => 3000, # 3000 kilometer minimum radius
LONGRID => 10, # longitude grid lines every 10 degrees
LATGRID => 10); # latitude grid lines every 10 degrees
#
## Plot points on the map
#
my ($lon, $lat) = getsomepoints();
$pl->map_plot ($lon, $lat, PLOTTYPE => 'POINTS', SYMBOL => 1, PROJECTION => 'AZEQDIST', CENTER => [-170, 70]);
#
## Plot lines on the map
#
my ($lon, $lat) = getsomemorepoints();
$pl->map_plot ($lon, $lat, PLOTTYPE => 'LINE',
PROJECTION => 'AZEQDIST',
CENTER => [-170, 70]);
For more information on GMT, see http://www.soest.hawaii.edu/gmt/
=head1 DESCRIPTION
This is the PDL/PLplot interface to the GMT map databases, allowing one to create
pleasing world maps in either of two projections:
1) linear (no projection)
2) Azimuthal Equidistant
The design is modular to allow addition of other projections.
Routines are also supplied to allow plotting of points and lines (supplied in degrees lon/lat)
on the maps with the correct projection supplied.
=head1 FUNCTIONS
=head2 worldmap
=for ref
Plot a world map using PLplot.
=for usage
Arguments: just a hash reference which can contain the following keywords:
PROJECTION : LINEAR (default) or AZEQDIST
For LINEAR projections:
MAPBOX: An array ref containing [WEST, EAST, SOUTH, NORTH] in degrees
-180 to 180, -90 to 90, ie: MAPBOX => [-180, 180, -90, 90]
For AZEQDIST projections:
CENTER : A list ref to the center point of the projection, in degrees, ie: [-170, 70]
RADIUS : A minimum radius in kilometers
For all projections:
RESOLUTION : The size of the map database used: "full", "high", "intermediate", "low" or "crude"
RIVER_DETAIL : A list reference to which rivers to plot:
1 = Permanent major rivers
2 = Additional major rivers
3 = Additional rivers
4 = Minor rivers
5 = Intermittent rivers - major
6 = Intermittent rivers - additional
7 = Intermittent rivers - minor
8 = Major canals
9 = Minor canals
10 = Irrigation canals
BOUNDARIES : A list reference to which boundaries to plot:
1 = National boundaries
2 = State boundaries within the Americas
3 = Marine boundaries
LONGRID : The grid spacing for longitude lines in degrees (undef = no lon grids)
LATGRID : The grid spacing for latitude lines in degrees (undef = no lat grids)
=head2 map_plot
=for ref
Plots lon/lat points or lines on an existing map with projection.
=for usage
map_plot ($lon, $lat, PROJECTION => ..., CENTER => [...,...]);
PROJECTION defaults to LINEAR. If AZEQDIST is specified, then the
CENTER lon/lat must be specified.
=head2 fetch
=for ref
Get lon and lat PDLs.
=for usage
Arguments:
A hash reference with these options available:
MAPBOX : An array ref containing [minlon, maxlon, minlat, maxlat] in degrees -180 to 180, -90 to 90
RESOLUTION : The size of the map database used: "full", "high", "intermediate", "low" or "crude"
RIVER_DETAIL : A list reference to which rivers to plot:
1 = Permanent major rivers
2 = Additional major rivers
3 = Additional rivers
4 = Minor rivers
5 = Intermittent rivers - major
6 = Intermittent rivers - additional
7 = Intermittent rivers - minor
8 = Major canals
9 = Minor canals
10 = Irrigation canals
BOUNDARIES : A list reference to which boundaries to plot:
1 = National boundaries
2 = State boundaries within the Americas
3 = Marine boundaries
Returns: ($lon, $lat) large 1-D PDLs
=for example
($lon, $lat) = PDL::Graphics::Map::fetch (MAPBOX => [-180, 180, -90, 90],
RESOLUTION => 'crude',
RIVER_DETAIL => [1,2,3,4]);
=head1 AUTHOR
Doug Hunt, dhunt\@ucar.edu.
=head1 SEE ALSO
perl(1), PDL(1), pscoast(l).
=cut
EOD
#-------------------------------------------------------------------------
# Perl portion of the interface (put by PP into the .pm file)
#-------------------------------------------------------------------------
pp_addpm (<<'EOPM');
use PDL::Primitive;
use PDL::Math;
use PDL::Core;
use PDL::Basic;
use PDL::Types;
use PDL::Slices;
use PDL::Graphics::PLplot;
use vars qw (%projection);
sub fetch {
my %parms = @_;
my @box = exists($parms{MAPBOX}) ? @{$parms{MAPBOX}} : (-180, 180, -90, 90);
die "bounding box must contain 4 edges: WESN in degrees (-180 -> 180, -90 -> 90)"
unless (@box == 4);
my $res = exists($parms{RESOLUTION}) ?
substr($parms{RESOLUTION}, 0, 1) : 'c'; # defaults to crude resolution
my @rivers = exists($parms{RIVER_DETAIL}) ? @{$parms{RIVER_DETAIL}} : ();
push (@rivers, (0) x 10); # note 10 river types
my @borders = exists($parms{BOUNDARIES}) ? @{$parms{BOUNDARIES}} : ();
push (@borders, (0) x 3); # note 3 boundary types
my $rlevels = pack ("i*", @rivers); # defaults to no rivers and canals
my $blevels = pack ("i*", @boundaries); # defaults to no national boundaries
my $drawc = 1; # defaults to 'draw coastlines'
my $lat = '';
my $lon = '';
pscoast($box[0], $box[1], $box[2], $box[3], $res, $rlevels, $blevels, $drawc, $lon, $lat);
my $size = length($lat)/8;
# Make a PDL of these data
my $latp = PDL->new; # Create piddle
$latp->set_datatype($PDL_D); # as a double array
$latp->setdims([$size]); # Set dimensions
${$latp->get_dataref} = $lat; # Assign the data
$latp->upd_data(); # Sync up everything - $cp is ready to be used.
my $lonp = PDL->new; # Create piddle
$lonp->set_datatype($PDL_D); # as a double array
$lonp->setdims([$size]); # Set dimensions
${$lonp->get_dataref} = $lon; # Assign the data
$lonp->upd_data(); # Sync up everything - $cp is ready to be used.
$lonp->inplace->setnantobad;
$latp->inplace->setnantobad;
return ($lonp, $latp);
}
# Convert lat/lon (degrees, -90 to 90, -180 to 180) to XY positions
# according to an azimuthal eqidistant projection
# (see http://mathworld.wolfram.com/StereographicProjection.html)
sub lonlat2azequi {
my $lon = shift; # PDL of one or more lons
my $lat = shift; # "" "" lats
my $lon0 = shift; # reference (center) longitude (scalar)
my $lat0 = shift; # reference (center) latitude
my $pi = 3.141592653589793238;
my ($a, $b) = (6378.1363, 6356.7516); # equatorial, polar Earth radii
my $r = ($a+$b)/2; # average Earth radius (use spherical approximate projections)
my $del = 1e-6;
my $m = ((abs($lat - $lat0) < $del) * (abs($lon - $lon0) < $del))->setbadtoval(0);
# subtract modulo 360 (acck!)
my $clon = $lon - $lon0;
while (any $clon > 180) { $clon = ($clon > 180) * ($clon-360) + ($clon <= 180) * $clon; }
while (any $clon < -180) { $clon = ($clon < -180) * ($clon+360) + ($clon >= -180) * $clon; }
# convert lat/lons to radians
$clon = $clon * ($pi/180);
$lat = $lat * ($pi/180);
$lon0 = pdl ($lon0 * ($pi/180))->dummy(0,$clon->nelem);
$lat0 = pdl ($lat0 * ($pi/180))->dummy(0,$clon->nelem);
my $c = acos ( sin($lat0)*sin($lat) + cos($lat0)*cos($lat)*cos($clon) );
my $k = $r * $c/sin($c);
my $x = $k * cos($lat)*sin($clon);
my $y = $k * (cos($lat0)*sin($lat) - sin($lat0)*cos($lat)*cos($clon) );
# set all points at origin to 0,0
(my $t = $x->where($m)) .= 0;
($t = $y->where($m)) .= 0;
# set any NaNs generated in the projection to the bad value
$x->inplace->setnantobad;
$y->inplace->setnantobad;
return ($x, $y);
}
# Convert XY positions to lat/lon (degrees, -90 to 90, -180 to 180)
# according to an azimuthal eqidistant projection
# (see http://mathworld.wolfram.com/StereographicProjection.html)
sub azequi2lonlat {
my $x = shift; # PDL of one or more x coordinates
my $y = shift; # "" "" y ""
my $lon0 = shift; # reference (center) longitude (scalar)
my $lat0 = shift; # reference (center) latitude
my $pi = 3.141592653589793238;
my ($a, $b) = (6378.1363, 6356.7516); # equatorial, polar Earth radii
my $r = ($a+$b)/2; # average Earth radius (use spherical approximate projections)
my $case = 1;
if ($lat0 == 90) { $case = 2; } elsif ($lat0 == -90) { $case = 3; }
$lon0 = pdl ($lon0 * ($pi/180))->dummy(0,$x->nelem);
$lat0 = pdl ($lat0 * ($pi/180))->dummy(0,$y->nelem);
$x /= $r;
$y /= $r;
my $c = sqrt($x**2 + $y**2);
my $lat = asin (cos($c)*sin($lat0) + ($y*sin($c)*cos($lat0))/$c);
my $lon;
if ($case == 1) {
$lon = $lon0 + atan2 ($x*sin($c), $c*cos($lat0)*cos($c) - $y*sin($lat0)*sin($c));
} elsif ($case == 2) {
$lon = $lon0 + atan2 ($x,-$y);
} else {
$lon = $lon0 + atan2 ($x, $y);
}
# convert to degrees
$lon = $lon * (180/$pi);
$lat = $lat * (180/$pi);
return ($lon, $lat);
}
# map of projection names to projection subroutines
%projection = (LINEAR => sub { return ($_[0], $_[1]); }, # lon/lat = x/y for LINEAR projection
AZEQDIST => \&lonlat2azequi);
# Draw points with projection
sub map_plot {
my $self = shift;
my $lon = shift;
my $lat = shift;
my %parms = @_;
my $proj = exists($parms{PROJECTION}) ? $parms{PROJECTION} : 'LINEAR'; # defaults to LINEAR
delete $parms{PROJECTION};
my @o = (0,0); # dummy projection center
if ($proj eq 'AZEQDIST') {
die "Must supply projection center point (CENTER = [lon, lat]) for AZEQDIST projection"
unless (exists($parms{CENTER}) && ((@o = @{$parms{CENTER}}) == 2));
delete $parms{CENTER};
}
# project lon/lat -> x/y
my ($x, $y) = &{$projection{$proj}}($lon, $lat, @o);
# plot
$self->xyplot ($x, $y, %parms);
}
# Draw a map of some section of the World in various projections and with various
# options. See POD doc above for details.
sub worldmap {
my $self = shift;
my %parms = @_;
# get rid of options not allowed in PLplot
my $proj = exists($parms{PROJECTION}) ? $parms{PROJECTION} : 'LINEAR'; # defaults to LINEAR
my @b = exists($parms{MAPBOX}) ? @{$parms{MAPBOX}} : (-180, 180, -90, 90); # bounding box in degrees
my $longrid = $parms{LONGRID};
my $latgrid = $parms{LATGRID};
my $resolution = $parms{RESOLUTION} || 'crude';
my $boundaries = exists ($parms{BOUNDARIES}) ? $parms{BOUNDARIES} : [1];
foreach (qw(PROJECTION MAPBOX LONGRID LATGRID RESOLUTION BOUNDARIES)) {
delete $parms{$_};
}
my @bxy = @b; # assume linear projection for now # bounding box in XY after projection
my @o = (0,0); # dummy projection center
my ($a, $b) = (6378.1363, 6356.7516); # equatorial and polar Earth radii
my $pi = 3.141592653589793238;
# compute map edges if only a center/radius is given (Azimuthal equidistant projections)
if ($proj eq 'AZEQDIST') {
die "Must supply projection center point (CENTER = [lon, lat]) for AZEQDIST projection"
unless (exists($parms{CENTER}) && ((@o = @{$parms{CENTER}}) == 2));
delete $parms{CENTER};
die "Must supply projection radius (RADIUS = val_in_km) for AZEQDIST projection"
unless (exists($parms{RADIUS}) && ((my $r = pdl($parms{RADIUS})) > 0));
delete $parms{RADIUS};
# determine edge points from center and radius
my ($lonb, $latb) = azequi2lonlat (append(-$r, $r), append(-$r, $r), @o);
# if box goes over the pole, set max lat to 90
if ($o[1]+($r/$a)*(180/$pi) > 80) {
$b[0] = -180;
$b[1] = 180;
$b[2] = $latb->min;
$b[3] = 90;
} elsif ($o[1]-($r/$a)*(180/$pi) < -80) { # set min lat to -90
$b[0] = -180;
$b[1] = 180;
$b[2] = -90;
$b[3] = $latb->max;
} else {
@b = ($lonb->list, $latb->list);
}
@bxy = (-$r, $r, -$r, $r);
} # end proj == AZEQDIST
#
## get coasts/rivers/borders from GMT database
#
my $minlon = $b[0] - 20 < -180 ? -180 : $b[0] - 20;
my $maxlon = $b[1] + 20 > 180 ? 180 : $b[1] + 20;
my $minlat = $b[2] - 10 < -90 ? -90 : $b[2] - 10;
my $maxlat = $b[3] + 10 > 90 ? 90 : $b[3] + 10;
my $box = [$minlon, $maxlon, $minlat, $maxlat];
my ($lonmap, $latmap) = fetch(%parms, BOUNDARIES => $boundaries, RESOLUTION => $resolution, MAPBOX => $box);
$lonmap->where($lonmap > 180) -= 360; # normalize to -180 to 180 longitudes D. Hunt 7/24/2008
delete $parms{RIVER_DETAIL}; # prevent error in xyplot later on. D. Hunt 7/24/2008
# project map lon/lat -> x/y
my ($xmap, $ymap) = &{$projection{$proj}}($lonmap, $latmap, @o);
# plot map
$self->xyplot ($xmap, $ymap,
BOX => [@bxy],
XBOX => 'BC', YBOX => 'BC',
JUST => 1,
LINEWIDTH => 1, COLOR => 'BLUE',
PLOTTYPE => 'LINE', %parms);
#
## now deal with lon/lat lines
#
my $n = 50; # the number of points plotted along the lat/lon lines
# compute longitude lines for map
if (defined($longrid)) {
my $lonlines = ((sequence(360/$longrid)*$longrid)-180);
my $nlonlines = $lonlines->nelem;
$lonlines = $lonlines->dummy(0,$n)->append(zeroes(1)/0)->clump(2); # put NaNs in to separate lines
# corresponding latitudes for longitude lines
my $lonlineslats = sequence($n)*(180/$n)-90;
$lonlineslats = $lonlineslats->dummy(0,$nlonlines)->xchg(0,1)->append(zeroes(1)/0)->clump(2);
# project to map
($lonlines, $lonlineslats) = &{$projection{$proj}}($lonlines, $lonlineslats, @o);
$lonlines->inplace->setnantobad;
$lonlineslats->inplace->setnantobad;
# get rid of lines off map
my $d = $bxy[3] - $bxy[2];
# mask of all lon values within box
my $m = (
($lonlineslats >= ($bxy[2] - 2*$d)) *
($lonlineslats <= ($bxy[3] + 2*$d))
);
$m->inplace->setbadtoval(1); # include bad values
$lonlines = $lonlines->where($m);
$lonlineslats = $lonlineslats->where($m);
# Plot
$self->xyplot ($lonlines, $lonlineslats,
XBOX => '', YBOX => '',
PLOTTYPE => 'LINE',
LINEWIDTH => 1, COLOR => 'BLACK', %parms) if (defined($lonlines));
#
## lon line labels
#
my $lablons = ((sequence(360/$longrid)*$longrid)-180);
my $lablats = pdl(($b[2] + $b[3])/2)->dummy(0,$lablons->nelem);
my $lablats1 = pdl(($b[2] + $b[3])/2.1)->dummy(0,$lablons->nelem);
# map projection
my ($x1, $y1) = &{$projection{$proj}}($lablons, $lablats, @o);
my ($x2, $y2) = &{$projection{$proj}}($lablons, $lablats1, @o);
my $dx = $x2-$x1;
my $dy = $y2-$y1;
# make sure the label is on the map
$m = (($x1 > $bxy[0]) * ($x1 < $bxy[1]) * ($y1 > $bxy[2]) * ($y1 < $bxy[3]));
# Plot longitude labels
$self->setparm(CHARSIZE => 0.85);
for (my $i=0;$i<$x1->nelem;$i++) {
next unless ($m->at($i));
$self->text (int($lablons->at($i)),
TEXTPOSITION => [$x1->at($i), $y1->at($i), $dx->at($i), $dy->at($i), 0.5], %parms);
}
}
# compute longitude lines for map
if (defined($latgrid)) {
# compute latitude lines for map
my $latlines = ((sequence(180/$latgrid)*$latgrid)-90);
$latlines = $latlines->where(abs($latlines) != 90); # get rid of poles
my $nlatlines = $latlines->nelem;
$latlines = $latlines->dummy(0,$n+1)->append(zeroes(1)/0)->clump(2); # put in NaNs in to separate lines
# corresponding longitudes for latitude lines
my $latlineslons = sequence($n+1)*(360/$n)-180;
$latlineslons = $latlineslons->dummy(0,$nlatlines)->xchg(0,1)->append(zeroes(1)/0)->clump(2);
# project to map
($latlineslons, $latlines) = &{$projection{$proj}}($latlineslons, $latlines, @o);
$latlines->inplace->setnantobad;
$latlineslons->inplace->setnantobad;
# This seems bogus and gets rid of BAD values which delimit separate sections.
# get rid of lines off map
my $d = $bxy[1] - $bxy[0];
my $m = ($latlineslons >= ($bxy[0] - 2*$d) * ($latlineslons <= ($bxy[1] + 2*$d))); # mask of all lon values within box
$m->inplace->setbadtoval(1); # include bad values
$latlines = $latlines->where($m);
$latlineslons = $latlineslons->where($m);
# Plot
$self->xyplot ($latlineslons, $latlines,
XBOX => '', YBOX => '',
PLOTTYPE => 'LINE',
LINEWIDTH => 1, COLOR => 'BLACK', %parms) if (defined($latlines));
#
## lat line labels
#
my $lablats = ((sequence(180/$latgrid)*$latgrid)-90);
my $lablons = pdl(($b[0] + $b[1])/2)->dummy(0,$lablats->nelem);
my $lablons1 = pdl(($b[0] + $b[1])/2.1)->dummy(0,$lablats->nelem);
# map projection
my ($x1, $y1) = &{$projection{$proj}}($lablons, $lablats, @o);
my ($x2, $y2) = &{$projection{$proj}}($lablons1,$lablats, @o);
my $dx = $x2-$x1;
my $dy = $y2-$y1;
# make sure the label is on the map
$m = (($x1 > $bxy[0]) * ($x1 < $bxy[1]) * ($y1 > $bxy[2]) * ($y1 < $bxy[3]));
# Plot latitude labels
$self->setparm(CHARSIZE => 0.85);
for (my $i=0;$i<$x1->nelem;$i++) {
next unless ($m->at($i));
$self->text (int($lablats->at($i)),
TEXTPOSITION => [$x1->at($i), $y1->at($i), $dx->at($i), $dy->at($i), 0.5], %parms);
}
}
}
EOPM
#-------------------------------------------------------------------------
# XS code for pscoast
#-------------------------------------------------------------------------
pp_addxs (<<'EOXS');
void
pscoast (west, east, south, north, res, rlevels, blevels, draw_coast, lon, lat)
double west
double east
double south
double north
char res
int *rlevels
int *blevels
int draw_coast
SV *lon
SV *lat
CODE:
{
pscoast (west, east, south, north, res, rlevels, blevels, draw_coast, lon, lat);
}
OUTPUT:
lon
lat
EOXS
pp_done();