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