GD::Map::Mercator
package GD::Map::Mercator;
use strict;
use warnings;
use GD;
use GD::Polyline;
our $VERSION = '1.03';
#
# keep map of regions to lat/long bounding boxes,
# so we can optimize dataset loading
#
our %regions = (
'africa-bdy' => [ -30.658889, -17.075556, 41.594722,53.089722],
'africa-cil' => [ -54.462778, -25.360556, 65.650278,77.588889],
'africa-riv' => [ -34.400278, -16.708611, 42.038056,55.2525],
'asia-bdy' => [ -9.126944, 19.001389, 54.476667,141.008889],
'asia-cil' => [ -54.753889, -190.351944, 81.851944,180.0],
'asia-riv' => [ -46.569722, -179.988056, 74.412222,180.0],
'europe-bdy' => [ 36.151389, -8.751389, 70.088889,31.586667],
'europe-cil' => [ 34.808889, -31.29, 71.165000,31.250278],
'europe-riv' => [ 36.566111, -21.791944, 69.999722,29.468889],
'namer-bdy' => [ 41.675556, -141.003056, 69.645556,-66.901111],
'namer-cil' => [ 24.538333, -168.132778, 83.623611,-12.155],
'namer-pby' => [ 29.706944, -139.047778, 68.900556,-57.1],
'namer-riv' => [ 25.768333, -166.053056, 74.032222,-54.470833],
'samer-bdy' => [ -55.556389, -117.1225, 32.718333,-51.682778],
'samer-cil' => [ -85.470278, -179.987778, 33.003333,179.976111],
'samer-riv' => [ -52.733333, -117.126111, 32.718333,-34.917222],
);
our %datasets = (
africa => [ 'bdy', 'cil', 'riv' ],
asia => [ 'bdy', 'cil', 'riv' ],
europe => [ 'bdy', 'cil', 'riv' ],
namer => [ 'bdy', 'cil', 'riv', 'pby' ],
samer => [ 'bdy', 'cil', 'riv' ],
);
our %colors = (
white => [255,255,255],
lgray => [191,191,191],
gray => [127,127,127],
dgray => [63,63,63],
black => [0,0,0],
lblue => [0,0,255],
blue => [0,0,191],
dblue => [0,0,127],
gold => [255,215,0],
lyellow => [255,255,0],
yellow => [191,191,0],
dyellow => [127,127,0],
lgreen => [0,255,0],
green => [0,191,0],
dgreen => [0,127,0],
lred => [255,0,0],
red => [191,0,0],
dred => [127,0,0],
lpurple => [255,0,255],
purple => [191,0,191],
dpurple => [127,0,127],
lorange => [255,183,0],
orange => [255,127,0],
pink => [255,183,193],
dpink => [255,105,180],
marine => [127,127,255],
cyan => [0,255,255],
lbrown => [210,180,140],
dbrown => [165,42,42],
transparent => [1,1,1]
);
my %valid_fmt = qw(png newFromPng gif newFromGif jpg newFromJpeg jpeg newFromJpeg);
sub new {
my ($class, %opts) = @_;
if ($opts{data_path}) {
die "data_path $opts{data_path} not found"
unless(-d $opts{data_path});
$opts{data_path} .= '/'
unless (substr($opts{data_path}, -1) eq '/');
}
if ($opts{basemap_path}) {
die "basemap_path $opts{basemap_path} not found."
unless(-d $opts{basemap_path});
$opts{basemap_path} .= '/'
unless (substr($opts{basemap_path}, -1) eq '/');
}
else {
$opts{basemap_path} = '';
}
my ($name, $fmt) = ('', 'png');
if ($opts{basemap_name}) {
$opts{basemap_name} .= '.png'
unless $opts{basemap_name}=~/\.(?:png|gif|jpg|jpeg)$/;
($name, $fmt) = ($opts{basemap_name}=~/^(.+)\.(png|gif|jpg|jpeg)$/);
}
$opts{silent} = 0 unless exists $opts{silent};
my $self = {
data_path => $opts{data_path},
basemap_path => $opts{basemap_path},
basemap_loc => $opts{basemap_path} . $name,
imgfmt => $fmt,
verbose => !$opts{silent},
thickness => $opts{thickness} || 2,
foreground => $opts{foreground} || 'black',
background => $opts{background} || 'white',
keeps => $opts{keep} || {},
omit => $opts{omit},
save_coords => $opts{save_coords},
};
bless $self, $class;
#
# if min/max info is provided, create a new basemap
# else attempt to load existing basemap
#
my $haspts = 0;
$haspts += defined($opts{$_}) ? 1 : 0
foreach (qw(min_lat min_long max_lat max_long));
die "Incomplete latitude/longitude datapoints provided."
if ($haspts > 0) && ($haspts < 4);
#
# attempt to load existing basemap; dies on any error
#
die "No basemap path or coordinates provided."
unless $haspts || ($self->{basemap_path} && $self->{basemap_loc});
return $self->_load_basemap()
unless $haspts;
#
# create a new basemap (maybe we should go ahead and check for a matching
# existing basemap ?)
#
die "No data_path specified."
unless $self->{data_path};
my ($minlat, $minlong, $maxlat, $maxlong, $width, $height) =
($opts{min_lat}, $opts{min_long},
$opts{max_lat}, $opts{max_long},
$opts{width} || 400, $opts{height} || 400);
my $mercator;
$self->{mercator} = $mercator = GD::Map::Mercator::Projector->new(
$minlat, $minlong, $maxlat, $maxlong, $width, $height, $self->{verbose});
my ($bg, $fg, $linew) = @$self{qw(background foreground thickness)};
$| = 1 if $self->{verbose};
#
# create empty image before loading data
# Note that image dimensions are adjusted by the projector
#
print "Creating GD image ($width x $height)\n"
if $self->{verbose};
($width, $height) = $mercator->dimensions();
my $img = $self->{image} = GD::Image->new($width, $height);
$self->{fg} = ref $fg ? $img->colorAllocate(@$fg) :
($fg=~/^#([0-9A-F]{2})([0-9A-F]{2})([0-9A-F]{2})$/i)
? $img->colorAllocate(hex($1), hex($2), hex($3))
: $img->colorAllocate(@{$colors{$fg}});
$self->{bg} = ref $bg ? $img->colorAllocate(@$bg) :
($bg=~/^#([0-9A-F]{2})([0-9A-F]{2})([0-9A-F]{2})$/i)
? $img->colorAllocate(hex($1), hex($2), hex($3))
: $img->colorAllocate(@{$colors{$bg}});
$img->filledRectangle(0,0,$width-1,$height-1,$self->{bg});
$img->setThickness($linew);
my @xy;
my $seg = [];
#
# check which of the datasets we need
#
my @regions = ();
my %omits = ();
map $omits{$_} = 1, @{$opts{omit}}
if $opts{omit};
foreach (keys %regions) {
if ((!$omits{substr($_, -3)}) &&
(
(($regions{$_}[0] <= $minlat) &&
($regions{$_}[2] >= $minlat)) ||
(($regions{$_}[0] <= $maxlat) &&
($regions{$_}[2] >= $maxlat)) ||
(($regions{$_}[0] >= $minlat) &&
($regions{$_}[2] <= $maxlat))
) &&
(
(($regions{$_}[1] <= $minlong) &&
($regions{$_}[3] >= $minlong)) ||
(($regions{$_}[1] <= $maxlong) &&
($regions{$_}[3] >= $maxlong)) ||
(($regions{$_}[1] >= $minlong) &&
($regions{$_}[3] <= $maxlong))
)) {
print "Using $_\n" if $self->{verbose};
push(@regions, $_);
}
}
# open each of the data files
my $oldfd = select(STDOUT);
$| = 1;
select $oldfd;
foreach (@regions) {
my $datapath = "$opts{data_path}$_.bin";
print "\nLoading datafile $datapath\n"
if $self->{verbose};
my $fd;
print "Skipping $datapath ($!)\n" and next
unless open $fd, $datapath;
binmode $fd;
$mercator->filter($fd, $self, $_, $self->{keep}{$_});
close $fd;
}
return $self->{basemap_loc}
? $self->save("$self->{basemap_loc}.$self->{imgfmt}")
: $self;
}
sub config {
return $_[0]->{mercator}->config();
}
sub extract {
my ($self, $minlat, $minlong, $maxlat, $maxlong, $scale) = @_;
my @coords = $self->{mercator}->config();
$@ = 'Specified region outside the bounds of this map.',
return undef
if ($minlat < $coords[0]) || ($maxlat > $coords[2]) ||
($minlong < $coords[1]) || ($maxlong > $coords[3]);
my ($xmin, $ymin, $xminmerc, $yminmerc) = $self->{mercator}->project($minlat, $minlong);
my ($xmax, $ymax, $xmaxmerc, $ymaxmerc) = $self->{mercator}->project($maxlat, $maxlong);
my ($width, $height) = (($xmax - $xmin + 1), ($ymax - $ymin + 1));
my $class = ref $self;
if ($scale) {
#
# just create a new map from scratch that scales the selected region
#
$width *= $scale;
$height *= $scale;
return ${class}->new(
data_path => $self->{data_path},
min_lat => $minlat,
min_long => $minlong,
max_lat => $maxlat,
max_long => $maxlong,
width => $width,
height => $height,
verbose => $self->{verbose},
thickness => $self->{thickness},
foreground => $self->{foreground},
background => $self->{background},
keep => $self->{keep},
omit => $self->{omit},
save_coords => $self->{save_coords},
);
}
my $newmap = {
data_path => $self->{data_path},
min_lat => $minlat,
min_long => $minlong,
max_lat => $maxlat,
max_long => $maxlong,
width => $width,
height => $height,
verbose => $self->{verbose},
thickness => $self->{thickness},
foreground => $self->{foreground},
background => $self->{background},
keep => $self->{keep},
omit => $self->{omit},
save_coords => $self->{save_coords},
};
$newmap->{mercator} = GD::Map::Mercator::Projector->new(
$minlat, $minlong, $maxlat, $maxlong, $width, $height, $self->{verbose});
($width, $height) = $newmap->{mercator}->dimensions();
my ($fg, $bg) = @$self{qw(foreground background)};
my $img = $newmap->{image} = GD::Image->new($width, $height);
$newmap->{fg} = ref $fg ? $img->colorAllocate(@$fg) :
($fg=~/^#([0-9A-F]{2})([0-9A-F]{2})([0-9A-F]{2})$/i)
? $img->colorAllocate(hex($1), hex($2), hex($3))
: $img->colorAllocate(@{$colors{$fg}});
$newmap->{bg} = ref $bg ? $img->colorAllocate(@$bg) :
($bg=~/^#([0-9A-F]{2})([0-9A-F]{2})([0-9A-F]{2})$/i)
? $img->colorAllocate(hex($1), hex($2), hex($3))
: $img->colorAllocate(@{$colors{$bg}});
$img->filledRectangle(1,1,$width,$height,$newmap->{bg});
$img->setThickness($newmap->{thickness});
$newmap->{image}->copy($self->{image}, 0, 0, $xmax, $ymax, $width, $height);
return bless $newmap, $class;
}
sub scale {
my ($self, $scale) = @_;
my $class = ref $self;
my @coords = $self->{mercator}->config();
return ${class}->new(
data_path => $self->{data_path},
min_lat => $coords[0],
min_long => $coords[1],
max_lat => $coords[2],
max_long => $coords[3],
width => $coords[-2] * $scale,
height => $coords[-1] * $scale,
silent => !$self->{verbose},
thickness => $self->{thickness},
foreground => $self->{foreground},
background => $self->{background},
keep => $self->{keep},
omit => $self->{omit},
save_coords => $self->{save_coords},
);
}
sub project {
my $self = shift;
return $self->{mercator}->project(@_);
}
sub translate {
my $self = shift;
return $self->{mercator}->translate(@_);
}
sub image { return $_[0]->{image}; }
sub save {
my ($self, $path) = @_;
my $fmt = (substr($path, -4, 1) eq '.') ? substr($path, -3)
: (substr($path, -5, 1) eq '.') ? substr($path, -4)
: undef;
my ($confpath, $imgpath, $coordpath) = ($path, $path, $path);
if ($fmt) {
$@ = "Invalid or unsupported image format $fmt",
return undef
unless $valid_fmt{$fmt} && $self->{image}->can($valid_fmt{$fmt});
$confpath=~s/$fmt$/conf/;
$coordpath=~s/$fmt$/coords/;
}
else {
$confpath .= '.conf';
$coordpath .= '.coords';
$imgpath .= '.png';
$fmt = 'png';
}
die "Cannot open $imgpath: $!" unless open OUTF, ">$imgpath";
binmode OUTF;
print "Writing $imgpath\n"
if $self->{verbose};
print OUTF ($fmt eq 'png') ? $self->{image}->png
: ($fmt eq 'gif') ? $self->{image}->gif
: $self->{image}->jpeg;
close OUTF;
die "Cannot open $confpath: $!"
unless open OUTF, ">$confpath";
#
# save as simple CSV, no D::D madness
#
print OUTF join(',', $self->{mercator}->config()), "\n";
close OUTF;
if ($self->{save_coords}) {
die "Cannot open $coordpath: $!"
unless open OUTF, ">$coordpath";
print OUTF $self->{_coords};
close OUTF;
}
return $self;
}
sub dimensions {
return $_[0]->{mercator}->dimensions();
}
############################
# PRIVATE FUNCTIONS
############################
#
# actually, this might be considered a pucliy overloaded
# function; its the callback from the Mercator::Projector
# object to render a segment
#
# NOTE: need to apply some compress here: reduce adjacent coords
# to a single pair of coords of equivalent line segment
#
sub add_segment {
my ($self, $seg, $dataset, $segno) = @_;
print "\rDrawing segment $segno "
if $self->{verbose};
if ($self->{save_coords}) {
my $coords = $self->{_coords} || '';
my $area = 0;
my $i = 0;
while ($i < @$seg) {
$i += 2
while ($i < @$seg) && (!defined $seg->[$i]);
last unless ($i < @$seg);
$coords .= "Dataset $dataset Segment $segno Area $area: ";
$coords .= join(',', $seg->[$i++], $seg->[$i++], '')
while ($i < @$seg) && (defined $seg->[$i]);
$coords .= "\n";
$area++;
}
$self->{_coords} = $coords;
}
my $fg = $self->{fg};
my $img = $self->{image};
my $polyline = GD::Polyline->new();
my ($lx, $ly) = ($seg->[0], $seg->[1]);
my $i = 2;
($lx, $ly) = ($seg->[$i++], $seg->[$i++])
unless defined $lx;
$polyline->addPt($lx, $ly);
my $ptcnt = 0;
while ($i < @$seg) {
my ($x, $y) = ($seg->[$i++], $seg->[$i++]);
$polyline->addPt($x, $y),
$ptcnt++,
($lx, $ly) = ($x, $y),
next
if defined $x;
#
# if an undef marker, draw current line
#
$img->polyline($polyline, $fg)
if ($ptcnt > 2);
$polyline = GD::Polyline->new();
($lx, $ly) = (undef, undef);
$ptcnt = 0;
next;
}
#
# draw any remaining line
#
$img->polyline($polyline, $fg)
if ($ptcnt > 1);
return 1;
}
sub _load_basemap {
my $self = shift;
my $pat = $valid_fmt{$self->{imgfmt}};
my $conf = "$self->{basemap_loc}.conf";
my $imgfile = "$self->{basemap_loc}.$self->{imgfmt}";
die "Unsupported image format $self->{imgfmt}; check your GD configuration."
unless GD::Image->can($pat);
die "$conf not found."
unless (-s $conf);
die "$imgfile not found."
unless (-s $imgfile);
die "Cannot open $conf: $!"
unless open F, $conf;
my $data = <F>;
close F;
chomp $data;
my @mapdata = split /,/, $data;
die "Invalid map data file"
unless (scalar @mapdata == 10);
#
# only install the lat/long and pixel coords, skip the mercator distances
#
$self->{mercator} = GD::Map::Mercator::Projector->new(
@mapdata[0..3], @mapdata[8..9], $self->{verbose});
my $fd;
die "Cannot open $imgfile: $!"
unless open $fd, $imgfile;
$self->{image} =
($self->{imgfmt} eq 'gif') ? GD::Image->newFromGif($fd)
: ($self->{imgfmt} eq 'png') ? GD::Image->newFromPng($fd)
: GD::Image->newFromJpeg($fd);
close $fd;
return $self;
}
1;
package GD::Map::Mercator::Projector;
use POSIX;
use Geo::Mercator;
use strict;
use warnings;
sub new {
my ($class, $minlat, $minlong, $maxlat, $maxlong, $width, $height, $verbose) = @_;
die "Bad min/max longitude"
if ($minlong < -180) || ($minlong > 180) ||
($maxlong < -180) || ($maxlong > 180) ||
(($maxlong < $minlong) &&
((($maxlong < 0) && ($minlong < 0)) ||
(($maxlong > 0) && ($minlong > 0))));
die "Bad min/max latitude"
if ($minlat > $maxlat) || ($minlat < -90) || ($minlat > 90) ||
($maxlat < -90) || ($maxlat > 90);
my ($xmin, $ymin) = mercate($minlat, $minlong);
my ($xmax, $ymax) = mercate($maxlat, $maxlong);
my $longadj = (($xmin > 0) && ($xmax < 0));
my $hscale = $width/($xmax - $xmin);
my $vscale = $height/($ymax - $ymin);
#
# adjust the image dimensions to match best scaling
#
my $scale = ($hscale < $vscale) ? $hscale : $vscale;
$height = _round(($ymax - $ymin) * $scale);
$width = _round(($xmax - $xmin) * $scale);
my ($minmerclong, $minmerclat) = mercate($minlat, $minlong);
my ($maxmerclong, $maxmerclat) = mercate($maxlat, $maxlong);
return bless {
minlat => $minlat,
minlong => $minlong,
maxlat => $maxlat,
maxlong => $maxlong,
minmerclat => $minmerclat,
minmerclong => $minmerclong,
maxmerclat => $maxmerclat,
maxmerclong => $maxmerclong,
xmin => $xmin,
ymin => $ymin,
xmax => $xmax,
ymax => $ymax,
scale => $scale,
longadj => $longadj,
width => $width,
height => $height,
verbose => $verbose}, $class;
}
sub config {
my $self = shift;
return (@$self{qw(minlat minlong maxlat maxlong
minmerclong minmerclat maxmerclong maxmerclat width height)});
}
sub dimensions { return ($_[0]->{width}, $_[0]->{height}); }
#
# return pixel coord from input lat/long
#
sub project {
my $self = shift;
#
# note: we assume the inputs are in a valid range, but not
# neccesarily inside our bounding box
#
my @result = ();
while (@_) {
my ($x, $y) = mercate(shift, shift);
push @result, (($y <= $self->{ymax}) &&
($y >= $self->{ymin}) &&
($x >= $self->{xmin}) &&
($x <= $self->{xmax}))
? (_round(($x - $self->{xmin}) * $self->{scale}),
# lat goes bottom to top, pixels top to bottom
_round($self->{height} - (($y - $self->{ymin}) * $self->{scale})))
: ();
}
return @result;
}
#
# return lat/long and mercator distances for input pixel coords
#
sub translate {
my $self = shift;
#
# upconvert to meters before demercating
#
my @result = ();
my ($xmin, $ymin, $scale, $longadj, $height) =
@$self{qw(xmin ymin scale longadj height)};
while (@_) {
my ($x, $y) = (shift, $height - shift);
$x /= $scale;
$x += $xmin;
$y /= $scale;
$y += $ymin;
push @result, demercate($x, $y), $y, $x;
}
return @result;
}
sub _round {
return (ceil($_[0]) - $_[0]) < ($_[0] - floor($_[0]))
? ceil($_[0])
: floor($_[0]);
}
sub filter {
my ($self, $fd, $container, $dataset, $keepers) = @_;
my ($xmin, $ymin, $xmax, $ymax, $scale, $longadj, $width, $height) =
@$self{qw(xmin ymin xmax ymax scale longadj width height)};
my ($n, $len, $record, $segno, $rank, $pts);
my @coords = ();
my @mercs = ();
my @keeppx = ();
#
# convert keep region coords to pixel coords
#
if ($keepers) {
my $i = 0;
my ($keepx, $keepy);
while ($i < @$keepers) {
($keepx, $keepy) = mercate($keepers->[$i++], $keepers->[$i++]);
push @keeppx,
(_round(($keepx - $xmin) * $scale),
_round($height - (($keepy - $ymin) * $scale)));
}
}
while ($n = read($fd, $len, 4)) {
$len = unpack('L', $len);
$n = read($fd, $record, $len);
($segno, $rank, $pts) = unpack('LLL', substr($record, 0, 12));
$pts *= 2;
#
# sometimes Perl doesn't believe me the 1st time!
#
@mercs = unpack("d$pts", substr($record, 12));
my $retry = 3;
$retry--,
@mercs = unpack("d$pts", substr($record, 12))
while ($retry && (@mercs != $pts));
die "No coords read!!!!" unless scalar @mercs == $pts;
print "\n*** Had to reload segment $segno ", 3 - $retry, " times!\n"
if $self->{verbose} && ($retry < 3);
my $inside = 0;
my ($lx, $ly) =
(_round(($mercs[0] - $xmin) * $scale),
_round($height - (($mercs[1] - $ymin) * $scale)));
my $i = 2;
while ($i < @mercs) {
my ($x, $y) =
(_round(($mercs[$i++] - $xmin) * $scale),
_round($height - (($mercs[$i++] - $ymin) * $scale)));
#
# scaling causes many pts to overlap, so skip them
# !!!NOTE need to optimze for pts on the same line segment
# (ie, no change in x, xor no change in y)
#
next
if ($x == $lx) && ($y == $ly);
if (($y <= $height) && ($y >= 0) && ($x >= 0) && ($x <= $width) &&
((!$keepers) || _keep(\@keeppx, $x, $y))) {
#
# if prior point outside, add it w/ a undef marker
# probably need to compute clipping intersection
#
push @coords, undef, undef, $lx, $ly
unless (($ly <= $height) && ($ly >= 0) && ($lx >= 0) && ($lx <= $width));
push @coords, $x, $y;
$inside++;
}
($lx, $ly) = ($x, $y);
}
if ($inside) {
$container->add_segment(\@coords, $dataset, $segno)
}
elsif ($self->{verbose}) {
print "\r*** Skipping segment $segno \r";
}
@mercs = ();
@coords = ();
}
print "\n" if $self->{verbose};
return $self;
}
sub _keep {
my ($keeppx, $x, $y) = @_;
my $i = 0;
$i += 4
while ($i < @$keeppx) &&
(($y > $keeppx->[$i+3]) || ($y < $keeppx->[1]) ||
($x < $keeppx->[$i]) || ($x > $keeppx->[$i+2]));
return ($i != @$keeppx);
}
1;
__END__