##############################
# Demo code to interpret the map data in earth.txt
#
# Craig DeForest, 17-Dec-2002
#
$| = 1;
print "Initializing...\n";
use PDL;
use PDL::NiceSlice;
# Snarf the file and build a separate piddle for each polyline
print "Interpreting map file...\n";
open(MAP,"<earth.txt");
$nelem = 0;
while(<MAP>) {
next if(m/^\#/ || m/^\s*$/);
s/^\s+//;
($x,$y,$z,$color) = split /\s+/;
if($color==0 and $#a > 1) {
push(@a,$a[0]);
my $c = $a[1]->[3];
print $c;
($a[-1])->[3] = $c;
$nelem++;
push(@mainlist,pdl(@a)); undef @a;
}
push(@a,[$x,$y,$z,$color ? $color-1:0]);
$nelem++;
}
print "Breaking up elements...\n";
$elements = zeroes(4,$nelem);
$pos = 0;
foreach $z(@mainlist) {
$n = $z->dim(1);
$elements->(0:2,$pos:$pos+$n-1) .= $z->(0:2);
$elements->((3),$pos:$pos+$n-2) .= $z->((3),1:-1);
$elements->((3),$pos+n-1) .= 0;
$pos += $n;
}
print "... $nelem vectors\n";
# Transform all the polylines into spherical coordinates.
use PDL::Transform;
$t = t_compose(t_scale(ones(3)*180/3.14159),t_spherical());
$lonlat = $t->apply( $elements->(pdl(2,0,1)) ) ->(1:2); # discard radius
# Clean up the pen values at the longitude singularity
print "cleaning up singularities...\n";
$p = $elements->((3));
$idx = which((abs($lonlat->((0),0:-2) - $lonlat->((0),1:-1))) > 355);
$p->($idx) .= 0 if($idx->nelem > 0);
# Plot map
use PDL::Graphics::PGPLOT::Window;
$w = pgwin(dev=>'/xs',size=>[10,10]);
$w->lines($lonlat,$p,{Title=>'Lat/Lon map of world coastlines',XTitle=>'East Longitude',YTitle=>'North Latitude',Axis=>2});
wfits($lonlat->glue(0,$p->dummy(0,1)),'earth_coast.vec.fits');
1;