1##############################
2# Demo code to interpret the map data in earth.txt
3#
4# Craig DeForest, 17-Dec-2002
5#
6$| = 1;
7print "Initializing...\n";
8use PDL;
9use PDL::NiceSlice;
10
11# Snarf the file and build a separate piddle for each polyline
12print "Interpreting map file...\n";
13
14open(MAP,"<earth.txt");
15$nelem = 0;
16while(<MAP>) {
17  next if(m/^\#/ || m/^\s*$/);
18  s/^\s+//;
19  ($x,$y,$z,$color) = split /\s+/;
20  if($color==0 and $#a > 1) {
21    push(@a,$a[0]);
22    my $c = $a[1]->[3];
23    print $c;
24    ($a[-1])->[3] = $c;
25    $nelem++;
26    push(@mainlist,pdl(@a)); undef @a;
27  }
28  push(@a,[$x,$y,$z,$color ? $color-1:0]);
29  $nelem++;
30}
31
32print "Breaking up elements...\n";
33$elements = zeroes(4,$nelem);
34$pos = 0;
35foreach $z(@mainlist) {
36  $n = $z->dim(1);
37  $elements->(0:2,$pos:$pos+$n-1) .= $z->(0:2);
38  $elements->((3),$pos:$pos+$n-2) .= $z->((3),1:-1);
39  $elements->((3),$pos+n-1) .= 0;
40  $pos += $n;
41}
42
43print "... $nelem vectors\n";
44
45# Transform all the polylines into spherical coordinates.
46use PDL::Transform;
47$t = t_compose(t_scale(ones(3)*180/3.14159),t_spherical());
48
49$lonlat = $t->apply( $elements->(pdl(2,0,1)) )  ->(1:2); # discard radius
50
51# Clean up the pen values at the longitude singularity
52print "cleaning up singularities...\n";
53$p = $elements->((3));
54$idx = which((abs($lonlat->((0),0:-2) - $lonlat->((0),1:-1))) > 355);
55$p->($idx) .= 0 if($idx->nelem > 0);
56
57# Plot map
58use PDL::Graphics::PGPLOT::Window;
59$w = pgwin(dev=>'/xs',size=>[10,10]);
60$w->lines($lonlat,$p,{Title=>'Lat/Lon map of world coastlines',XTitle=>'East Longitude',YTitle=>'North Latitude',Axis=>2});
61
62wfits($lonlat->glue(0,$p->dummy(0,1)),'earth_coast.vec.fits');
631;
64