plot-latlong/plot-latlong

   1 #!/usr/bin/perl -w
   2
   3 ##############################################################################
   4 ## Plots user supplied lat/longs on a user chosen map.
   5 ##
   6 ## plot-latlong reads lat/long values from stdin and writes a PNG to
   7 ## stdout.  The input should contain one pair of lat/long values per line,
   8 ## with the values separated by whitespace.  Blank lines and lines starting
   9 ## with the pound character ('#') are ignored.
  10 ##
  11 ## Command line options are
  12 ##
  13 ##  -m <map-name> specifies the name of a map to use (default: the first
  14 ##     map listed in .mapinfo: 'World'); see .mapinfo for the valid map names;
  15 ##     see the file CONFIG for a description of the .mapinfo format
  16 ##
  17 ##  -s <point-size> specifies the size of the points to draw (default: 1);
  18 ##     points are filled squares, and the size is the width in pixels
  19 ##
  20 ##  -c causes the pixel coordinates of each lat/long to be printed to stderr;
  21 ##     the coordinates (0, 0) are at the upper left corner, and values increase
  22 ##     to the right and down
  23 ##
  24 ##  -i <map-info-file> specifies an alternate location for .mapinfo (other
  25 ##     than in the current directory or $HOME).
  26 ##
  27 ##----------------------------------------------------------------------------
  28 ##
  29 ## The code for handling the Alber/Lambert map projection is derived from
  30 ## GTrace v1.0.0beta (http://www.caida.org/tools/visualization/gtrace),
  31 ## which was written by Ram Periakaruppan.  The included set of maps are also
  32 ## derived from the GTrace distribution.  GTrace redistributed these maps
  33 ## with the permission of VisualRoute (http://www.visualroute.com),
  34 ## the original source of the maps.
  35 ##
  36 ##----------------------------------------------------------------------------
  37 ##
  38 ## Copyright (C) 2003,2004,2005 The Regents of the University of California.
  39 ## All Rights Reserved.
  40 ##
  41 ## Permission to use, copy, modify and distribute any part of this
  42 ## plot-latlong software package for educational, research and non-profit
  43 ## purposes, without fee, and without a written agreement is hereby
  44 ## granted, provided that the above copyright notice, this paragraph
  45 ## and the following paragraphs appear in all copies.
  46 ##
  47 ## Those desiring to incorporate this into commercial products or use
  48 ## for commercial purposes should contact the Technology Transfer
  49 ## Office, University of California, San Diego, 9500 Gilman Drive, La
  50 ## Jolla, CA 92093-0910, Ph: (858) 534-5815, FAX: (858) 534-7345.
  51 ##
  52 ## IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY
  53 ## PARTY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL
  54 ## DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS
  55 ## SOFTWARE, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN ADVISED OF
  56 ## THE POSSIBILITY OF SUCH DAMAGE.
  57 ##
  58 ## THE SOFTWARE PROVIDED HEREIN IS ON AN "AS IS" BASIS, AND THE
  59 ## UNIVERSITY OF CALIFORNIA HAS NO OBLIGATION TO PROVIDE MAINTENANCE,
  60 ## SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. THE UNIVERSITY
  61 ## OF CALIFORNIA MAKES NO REPRESENTATIONS AND EXTENDS NO WARRANTIES
  62 ## OF ANY KIND, EITHER IMPLIED OR EXPRESS, INCLUDING, BUT NOT LIMITED
  63 ## TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A
  64 ## PARTICULAR PURPOSE, OR THAT THE USE OF THE SOFTWARE WILL NOT INFRINGE
  65 ## ANY PATENT, TRADEMARK OR OTHER RIGHTS.
  66 ##
  67 ## The plot-latlong software is developed by the plot-latlong Team at the
  68 ## University of California, San Diego under the Cooperative Association
  69 ## for Internet Data Analysis (CAIDA) Program.  Support for this work is
  70 ## provided by the National Communications System (NCS) via NSF grant
  71 ## ANI-0221172, entitled "Routing and Peering Analysis for Enhancing
  72 ## Internet Performance and Security."
  73 ##
  74 ##############################################################################
  75
  76 use strict;
  77 use File::Basename;
  78 use GD;
  79 use Getopt::Std;
  80 use vars qw($opt_m $opt_s $opt_c $opt_i);
  81
  82 sub usage
  83 {
  84   die "usage: cat datapoints | plot-latlong [-m <map-name>] [-s <point-size>] [-c] [-i <map-info-file>] >output.png\n";
  85 }
  86
  87 if (not getopts('m:s:ci:'))
  88 {
  89   usage();
  90 }
  91
  92 my $DEBUG = 0;
  93 my $PI = 3.141592654;
  94
  95 my $point_size = $opt_s || 1;
  96
  97 my $first_map;
  98 my $map_directory;  # directory containing map images
  99 my %configuration;
 100 my @mapinfo_locations = (".mapinfo", "$ENV{HOME}/.mapinfo");
 101 if ($opt_i) {
 102     unshift @mapinfo_locations, $opt_i;
 103 }
 104 load_configuration(@mapinfo_locations);
 105
 106 my $selected_map = $opt_m || $first_map;
 107
 108 # lat/long of the upper-left and lower-right corners of the map image
 109 my $map_top_lat;
 110 my $map_top_long;
 111 my $map_bottom_lat;
 112 my $map_bottom_long;
 113
 114 # parameters for linear projection: pixels per degree
 115 my $map_lat_scale;
 116 my $map_long_scale;
 117
 118 # parameters for Alber projection (supplied by user):
 119 # (in degrees)
 120 my $R;              # radius of sphere
 121 #my $phi_1;          # standard parallel
 122 #my $phi_2;          # standard parallel
 123 #my $phi_0;          # origin latitude
 124 #my $lambda_0;       # origin longitude
 125 my $false_easting;  # the false easting amount
 126 my $false_northing; # the false northing amount
 127
 128 # values for Alber projection calculated from user parameters:
 129 # (in radians)
 130 my $PHI_1;
 131 my $PHI_2;
 132 my $PHI_0;
 133 my $LAMBDA_0;
 134 my $N;
 135 my $C;
 136
 137 my $conversion_fn;  # pointer to function for converting lat/long to (x, y)
 138
 139 my $map = load_map($selected_map);
 140
 141 #############################################################################
 142
 143 my $green = $map->colorAllocate(64, 192, 64);
 144 my $red = $map->colorAllocate(220, 64, 64);
 145
 146 while (<>)
 147 {
 148   chomp;
 149   next if /^\s*$/;
 150   next if /^\s*#/;
 151
 152   # REQUIRE: -90 <= $lat <= 90
 153   my ($lat, $long) = split /\s+/;
 154   my $adjusted_long = $long;
 155   if ($long < $map_top_long) {
 156       $adjusted_long += 360;
 157   } elsif ($long > $map_bottom_long) {
 158       $adjusted_long -= 360;
 159   }
 160   my ($x, $y) = &$conversion_fn($lat, $adjusted_long);
 161
 162   print STDERR "$lat $long $x $y\n" if $DEBUG || $opt_c;
 163
 164   if ($point_size == 1)
 165   {
 166     $map->setPixel($x, $y, $red);
 167   }
 168   else
 169   {
 170     my $half = int($point_size / 2);
 171     my $left_x = $x - $half;
 172     my $top_y = $y - $half;
 173     $map->filledRectangle($left_x, $top_y,
 174                           $left_x + $point_size - 1, $top_y + $point_size - 1,
 175                           $red);
 176   }
 177 }
 178
 179 print $map->png if not $DEBUG;
 180
 181 #############################################################################
 182 #############################################################################
 183
 184 sub convert_latlong_to_xy_linear
 185 {
 186   my ($lat, $long) = @_;
 187
 188   my $lat_rel = $map_top_lat - $lat;
 189   my $long_rel = $long - $map_top_long;
 190
 191   return ($long_rel * $map_long_scale, $lat_rel * $map_lat_scale);
 192 }
 193
 194 #############################################################################
 195
 196
 197 # See the publication "Map Projections Used by the U.S. Geological
 198 # Survey Bulletin 1532" for details about this projection.
 199 #
 200 # However, the present coder does not know exactly which projection this
 201 # function corresponds to in the USGS Bulletin.
 202
 203 sub round {
 204     my ($number) = @_;
 205     return int($number + .5 * ($number <=> 0));
 206 }
 207
 208 sub convert_latlong_to_xy_alber
 209 {
 210   my ($lat, $long) = @_;
 211
 212   my $phi = ($lat * $PI) / 180.0;
 213   my $lambda = ($long * $PI) / 180.0;
 214
 215   my $p   = ($R * sqrt($C - 2.0 * $N * sin($phi))) / $N;
 216   my $p_0 = ($R * sqrt($C - 2.0 * $N * sin($PHI_0))) / $N;
 217   my $theta = $N * ($lambda - $LAMBDA_0);
 218
 219   my $x = $false_easting + round($p * sin($theta));
 220   my $y = $false_northing - round($p_0 - $p * cos($theta));
 221
 222   return ($x, $y);
 223 }
 224
 225 ############################################################################
 226
 227 sub load_map
 228 {
 229   my ($name) = @_;
 230
 231   if (not exists $configuration{"MAP $name"})
 232   {
 233     die "ERROR: Map '$name' not found in map configuration file.\n";
 234   }
 235
 236   print("OPEN: $name ", $configuration{"MAP $name"}, "\n") if $DEBUG;
 237
 238   my $filename;
 239   ($filename, $map_top_lat, $map_top_long, $map_bottom_lat, $map_bottom_long)
 240     = split(/\s+/, $configuration{"MAP $name"});
 241
 242   my $path = "$map_directory/$filename";
 243   my $retval = new GD::Image($path)
 244     or die "ERROR: Couldn't open map image '$path': $!\n";
 245
 246   my ($map_width, $map_height) = $retval->getBounds();
 247
 248   print("DIM: $map_width $map_height\n") if $DEBUG;
 249
 250   $map_lat_scale = $map_height / ($map_top_lat - $map_bottom_lat);
 251   $map_long_scale = $map_width / ($map_bottom_long - $map_top_long);
 252
 253   print("SCALE: $map_lat_scale $map_long_scale\n") if $DEBUG;
 254
 255   # -- projections -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 256
 257   if (not exists $configuration{"PROJECTION $name"})
 258   {
 259     $conversion_fn = \&convert_latlong_to_xy_linear;
 260   }
 261   else
 262   {
 263     my @F = split(/\s+/, $configuration{"PROJECTION $name"});
 264     if ($F[0] eq "ALBER")
 265     {
 266       ($R, $false_easting, $false_northing) = @F[1,6,7];
 267       my ($phi_1, $phi_2, $phi_0, $lambda_0) = @F[2..5];
 268
 269       $PHI_1 = ($phi_1 * $PI) / 180.0;
 270       $PHI_2 = ($phi_2 * $PI) / 180.0;
 271       $PHI_0 = ($phi_0 * $PI) / 180.0;
 272       $LAMBDA_0 = ($lambda_0 * $PI) / 180.0;
 273
 274       $N = (sin($PHI_1) + sin($PHI_2)) / 2.0;
 275       $C = cos($PHI_1) ** 2 + 2.0 * $N * sin($PHI_1);
 276
 277       $conversion_fn = \&convert_latlong_to_xy_alber;
 278     }
 279     else
 280     {
 281       die "INTERNAL ERROR: unknown projection '$F[0]'\n";
 282     }
 283   }
 284
 285   return $retval;
 286 }
 287
 288 #############################################################################
 289
 290 sub load_configuration
 291 {
 292   foreach my $filename (@_)
 293   {
 294     if (-f $filename)
 295     {
 296       $map_directory = dirname($filename) . "/.mapimages";
 297       load_configuration_aux($filename);
 298       return;
 299     }
 300   }
 301
 302   die "ERROR: Couldn't find map configuration file; looked for: @_\n";
 303 }
 304
 305 sub load_configuration_aux
 306 {
 307   my ($filename) = @_;
 308
 309   open CONFIG, "$filename"
 310     or die "ERROR: Couldn't open map configuration file '$filename': $!\n";
 311   while (<CONFIG>)
 312   {
 313     chomp;
 314     next if /^\s*$/;
 315     next if /^\s*#/;
 316
 317     my @F = split /\s+/;
 318     if ($F[0] eq "MAP")
 319     {
 320       if ($#F == 6)
 321       {
 322         my ($top_lat, $top_long, $bot_lat, $bot_long) = @F[3..6];
 323         if (!check_coordinate($top_lat)
 324             || !check_coordinate($top_long)
 325             || !check_coordinate($bot_lat)
 326             || !check_coordinate($bot_long))
 327         {
 328           die "ERROR: Line $.: boundary coordinates are malformed in map configuration file.\n";
 329         }
 330
 331         if ($top_lat < $bot_lat || $top_long > $bot_long)
 332         {
 333           die "ERROR: Line $.: boundary coordinates have wrong relations in map configuration file.\n";
 334         }
 335
 336         $first_map = $F[1] if not defined $first_map;
 337         $configuration{"MAP $F[1]"} = join(" ", @F[2..6]);
 338         next;
 339       }
 340     }
 341     elsif ($F[0] eq "PROJECTION")
 342     {
 343       if ($#F >= 2)
 344       {
 345         if ($F[2] eq "ALBER")
 346         {
 347           if ($#F == 9)
 348           {
 349             foreach my $x (@F[3..9])
 350             {
 351               if (not check_number($x))
 352               {
 353                 die "ERROR: Line $.: projection parameters are malformed in map configuration file.\n";
 354               }
 355             }
 356
 357             $configuration{"PROJECTION $F[1]"} = join(" ", @F[2..9]);
 358             next;
 359           }
 360         }
 361         else
 362         {
 363           die "ERROR: Line $.: unknown map projection '$F[2]' in map configuration file.\n";
 364         }
 365       }
 366     }
 367     else
 368     {
 369       die "ERROR: Line $.: unknown key '$F[0]' in map configuration file.\n";
 370     }
 371
 372     die "ERROR: Line $. of map configuration file is malformed.\n";
 373   }
 374   close CONFIG;
 375 }
 376
 377 #############################################################################
 378
 379 sub check_coordinate
 380 {
 381   my ($x) = @_;
 382
 383   return $x =~ /^(\-?)\d+(\.\d+)?$/;
 384 }
 385
 386
 387 sub check_number
 388 {
 389   return check_coordinate(@_);
 390 }