The Perl Toolchain Summit needs more sponsors. If your company depends on Perl, please support this very important event.
use Config;

pp_add_exported('', 'worldmap map_line map_points'); 

pp_addpm({At => Top}, <<'EOD');
=head1 NAME 

PDL::Graphics::PGPLOT::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 (cvs) version of PDL.
       Set WITH_BADVAL => 1 in perldl.conf).

=head1 SYNOPSIS

  use PDL;
  use PDL::Graphics::Map;
  use PDL::Graphics::PGPLOT;

  #
  ## plot just the coastlines in a linear projection
  #

  dev "pdlmap1.ps/CPS";
  worldmap ({WEST  => -180, EAST  => 180, SOUTH => -90, NORTH => 90});

  # or, alternatively, in two steps:
  dev "pdlmap2.ps/CPS";
  my ($lon, $lat) = PDL::Graphics::PGPLOT::Map::fetch (
                                     {BOX => [-180, 180, -90, 90],
                                      SEPARATOR    => -999});
  # Plot the entire set of coastline polylines 
  line $lon, $lat, {MISSING => -999};

  #
  ## Plot coastlines, lat/lon grids and labels in an Azimuthal Equidistant projection
  #

  dev "pdlmap3.ps/CPS";
  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();
  map_points ($lon, $lat, {PROJECTION => 'AZEQDIST', CENTER  => [-170, 70]});

  #
  ## Plot lines on the map
  #

  my ($lon, $lat) = getsomemorepoints();
  map_line ($lon, $lat, {PROJECTION => 'AZEQDIST', CENTER  => [-170, 70], MISSING => -999});
  
For more information on GMT, see http://www.soest.hawaii.edu/gmt/ 

=head1 DESCRIPTION

This is the PDL 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 

using PGPLOT.

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

=for usage

Arguments: just a hash reference which can contain the following keywords:

  PROJECTION  : LINEAR (default) or AZEQDIST

For LINEAR projections:

  WEST  : the Western border of the map (degrees west -180 - 180)  
  EAST  : the Eastern border of the map (degrees west -180 - 180)
  SOUTH : the Southern border of the map (degrees north -90 - 90)
  NORTH : the Northern border of the map (degrees north -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                                                           

  COASTS : A boolean value:  plot coasts = true, don't = false

  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_points

=for ref

Plots lon/lat points on an existing map with projection.

=for usage

  map_points ($lon, $lat, {PROJECTION => ..., CENTER => [...,...]});

  PROJECTION defaults to LINEAR.  If AZEQDIST is specified, then the 
  CENTER lon/lat must be specified.

=head2 map_line

=for ref

Plots lon/lat lines on an existing map with projection.

=for usage

  map_line ($lon, $lat, {PROJECTION => ..., CENTER => [...,...], MISSING => ...});

  PROJECTION defaults to LINEAR.  If AZEQDIST is specified, then the 
  CENTER lon/lat must be specified.  To plot more than one line
  segment, specify MISSING to be a separator value.

=head2 fetch

=for ref

Get lon and lat PDLs.

=for usage      	

Arguments:  
  A hash reference with these options available:
  BOX   : 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                                                           

  COASTS : A boolean value:  plot coasts = true, don't = false
  SEPARATOR : A numeric value:  The value to place between each separate line segment.

Returns:  ($lon, $lat) large 1-D PDLs

=for example
  ($lon, $lat) = PDL::Graphics::Map::fetch (
                                     {BOX => [-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

$VERSION = '0.09';

#-------------------------------------------------------------------------
# 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::PGPLOT;
use PGPLOT;
use vars qw (%projection);

sub fetch {
  my $parms   = shift;

  my @box = exists($$parms{BOX}) ? @{$$parms{BOX}} : (-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 $separator = exists($$parms{SEPARATOR}) ? $$parms{SEPARATOR} : -999;

  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 = $lonp->badmask($separator);
  $latp = $latp->badmask($separator);
  
  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;

  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_points {
  my $lon = shift;
  my $lat = shift;
  my $parms = shift;

  my $sym;
  unless (ref($parms) =~ /HASH/) {
    $sym   = $parms;
    $parms = shift;
  }

  my $proj = exists($$parms{PROJECTION}) ? $$parms{PROJECTION} : 'LINEAR';  #  defaults to LINEAR
  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));
  }

  # project lon/lat -> x/y
  my ($x, $y) = &{$projection{$proj}}($lon, $lat, @o);

  # plot
  points $x, $y, $sym, $parms;
}


# Draw lines with projection
sub map_line {
  my $lon = shift;
  my $lat = shift;
  my $parms = shift;

  my $proj = exists($$parms{PROJECTION}) ? $$parms{PROJECTION} : 'LINEAR';  #  defaults to LINEAR
  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));
  }

  # project lon/lat -> x/y
  my ($x, $y) = &{$projection{$proj}}($lon, $lat, @o);

  # plot
  line $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 $parms = shift;

  my $proj = exists($$parms{PROJECTION}) ? $$parms{PROJECTION} : 'LINEAR';  #  defaults to LINEAR
  my @b    = exists($$parms{BOX}) ? @{$$parms{BOX}} : (-180, 180, -90, 90); #  bounding box in degrees
  my @bxy  = @b;      # assume linear projection for now                    #  bounding box in XY after projection
  my @o    = (0,0);   # dummy projection center

  my $bad     = -99999;
  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));

    die "Must supply projection radius (RADIUS = val_in_km) for AZEQDIST projection"
     unless (exists($$parms{RADIUS}) && ((my $r = $$parms{RADIUS}) > 0));

    # 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


  #
  ## set up PGPLOT
  #

  my $win = PDL::Graphics::PGPLOT::get_current_window();
  warn "Unable to get current window handle.  Please upgrade your PDL installation" unless defined($win);
  if (!defined($win) || !$win->{Hold}) {  # set up plotting environment if not on hold
    pgscr 0, 1,1,1;    # set background color = white
    pgscr 3, 0,0,0;    # set axis color = black
    pgenv @bxy, 1, -2; # set bounds, specify equal scale on X/Y axes
    pgsci (3);         # set to black
    pgbox ('BC', 0, 0, 'BC', 0, 0);  # plot minimal border
    $win->hold();
  }

  return if ($$parms{SetupMap});  # don't do anything more if setup only requested

  #
  ## 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, BOX => $box, SEPARATOR  => $bad});

  # set separator to the bad value, so projection won't operate on it.
  $lonmap->inplace->setvaltobad($bad);
  $latmap->inplace->setvaltobad($bad);

  # project map lon/lat -> x/y
  my ($xmap, $ymap) = &{$projection{$proj}}($lonmap, $latmap, @o);

  # set any NaNs generated in the projection to the bad value
  $xmap->inplace->setnantobad;
  $ymap->inplace->setnantobad;
  $xmap = $xmap->setbadtoval($bad);
  $ymap = $ymap->setbadtoval($bad);

  # plot
  line $xmap, $ymap, {MISSING => $bad, LINEWIDTH => 1, COLOR => 4};

  # 
  ## now deal with lon/lat lines
  #

  my $n = 50;       # the number of points plotted along the lat/lon lines 

  # compute longitude lines for map
  my $loninc = $$parms{LONGRID};  
  if (defined($loninc)) {
    my $lonlines =  ((sequence(360/$loninc)*$loninc)-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)); 
    $lonlines     = $lonlines->where($m);
    $lonlineslats = $lonlineslats->where($m);

    $lonlines->inplace->badmask($bad);
    $lonlineslats->inplace->badmask($bad);

    # Plot
    line $lonlines, $lonlineslats, {MISSING => $bad, LINEWIDTH => 1, COLOR => 3} if (defined($lonlines));

    #
    ## lon line labels
    #

    my $lablons  = ((sequence(360/$loninc)*$loninc)-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 $angles = atan2(($y2-$y1),($x2-$x1)) * (180/$pi);
    $angles = ($angles >  90) * ($angles - 180) + ($angles <=  90) * $angles;
    $angles = ($angles < -90) * ($angles + 180) + ($angles >= -90) * $angles;

    # make sure the label is on the map
    $m = ($x1 > $bxy[0] & $x1 < $bxy[1] & $y1 > $bxy[2] & $y1 < $bxy[3]);

    # Plot longitude labels
    pgsch(0.85);  # small characters
    for (my $i=0;$i<$x1->nelem;$i++) {
      next unless ($m->at($i));
      pgptxt ($x1->at($i), $y1->at($i), $angles->at($i), 0.5, int($lablons->at($i)));
    }

  }


  # compute longitude lines for map
  my $latinc = $$parms{LATGRID};  
  if (defined($latinc)) {

    # compute latitude lines for map
    my $latlines     =  ((sequence(180/$latinc)*$latinc)-90);
    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;

    # 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 
    $latlines     = $latlines->where($m);
    $latlineslons = $latlineslons->where($m);

    $latlines->inplace->badmask($bad);
    $latlineslons->inplace->badmask($bad);

    # Plot
    line $latlineslons, $latlines, {MISSING => $bad, LINEWIDTH => 1, COLOR => 3} if (defined($latlines));

    #
    ## lat line labels
    #

    my $lablats  = ((sequence(180/$latinc)*$latinc)-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 $angles = atan2(($y2-$y1),($x2-$x1)) * (180/$pi);
    $angles = ($angles >  90) * ($angles - 180) + ($angles <=  90) * $angles;
    $angles = ($angles < -90) * ($angles + 180) + ($angles >= -90) * $angles;

    # make sure the label is on the map
    $m = ($x1 > $bxy[0] & $x1 < $bxy[1] & $y1 > $bxy[2] & $y1 < $bxy[3]);

    # Plot latitude labels
    pgsch(0.85);  # small characters
    for (my $i=0;$i<$x1->nelem;$i++) {
      next unless ($m->at($i));
      pgptxt ($x1->at($i), $y1->at($i), $angles->at($i), 0.5, int($lablats->at($i)));
    }

  } 

}            

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();