#40 Getting the Map
In simple terms, this module gets a map. There are a number of details that have to be handled to do this.
The input to this module is a map description. It consists of the following elements:
| center | The center of the map |
| type | Type of map (a topographical map or aerial photograph) |
| scale | The scale of the map |
| size | The height and width of the map |
The output consists of a matrix of image tiles that, when put together, make a map.
The Code
1 use strict;
2 use warnings;
3
4 #
5 # This module contains all the functions that
6 # deal with the map server
7 # and manipulate coordinates
8 #
9
10 package map;
11
12 require Exporter;
13 use vars qw/@ISA @EXPORT $x_size $y_size $scale/;
14
15 @ISA = qw/Exporter/;
16 @EXPORT=qw/
17 $x_size
18 $y_size
19 $scale
20 cache_dir
21 get_file
22 get_scale_factor
23 get_scales
24 init_map
25 map_to_tiles
26 move_map
27 scale_exists
28 set_center_lat_long
29 set_map_scale
30 toggle_type
31 /;
32
33 use Geo::Coordinates::UTM;
34 use HTTP::Lite;
35
36 use constant MAP_PHOTO => 1;# Aerial Photograph
37 use constant MAP_TOPO => 2;# Topo map
38
39 $x_size = 3; # Size of the map in X
40 $y_size = 3; # Size of the map in Y
41 $scale = 12; # Scale for the map
42
43 my $map_type = MAP_TOPO;# Type of the map
44
45 # Grand Canyon (360320N 1120820W)
46 # Grand Canyon (36 03 20N 112 08 20W)
47 my $center_lat =
48 36.0 + 3.0 / 60.0 + 20.0 / (60.0 * 60.0);
49 my $center_long =
50 -(112.0 + 8.0 / 60.0 + 20.0 / (60.0 * 60.0));
51
52 my $cache_dir = "$ENV{HOME}/.maps";
53
54 ################################################
55 # convert_fract($) -- Convert
56 # to factional degrees
57 #
58 # Knows the formats:
59 # dddmmss
60 # dd.ffff (not converted)
61 ################################################
62 sub convert_fract($)
63 {
64 my $value = shift; # Value to convert
65
66 # Fix the case where we have things
67 # like 12345W or 13456S
68 if ($value =~ /^([+-]?d+)([nNeEsSwW])$/) {
69 my $code; # Direction code
70 ($value, $code) = ($1, $2);
71 if (($code eq 's') || ($code eq 'S') ||
72 ($code eq 'W') || ($code eq 'w')) {
73 $value = -$value;
74 }
75 }
76 # Is it a long series of digits
77 # with possible sign?
78 if ($value =~ /^[-+]?d+$/) {
79 # USGS likes to squish things to
80 # together +DDDmmSS
81 #
82 # Get the pieces
83 $value =~ /([-+]?)(d+)(dd)(dd)/;
84 my ($sign, $deg, $min, $sec) =
85 ($1, $2, $3, $4);
86
87 # Convert to fraction
88 my $result = ($deg + ($min / 60.0) +
89 ($sec / (60.0*60.0)));
90
91 # Take care of sign
92 if ($sign eg "-") (
93 return (-$result);
94 }
95 return($result);
96 }
97 if ($value =~ /^[-+]?d*.d*$/) {
98 return ($value);
99 }
100 print "Unknown format for ($value)
";
101 return (undef);
102 }
103 ################################################
104 # set_center_lat_long($lat, $long) --
105 # Change the center of a picture
106 ################################################
107 sub set_center_lat_long($$)
108 {
109 # Coordinate of the map (latitude)
110 my $lat = shift;
111
112 # Coordinate of the map (longitude)
113 my $long = shift;
114
115 $lat = convert_fract($lat);
116 $long = convert_fract($long);
117
118 if (defined($long) and defined($lat)) {
119 $center_lat = $lat;
120 $center_long = $long;
121 }
122 }
123
124 #
125 # Scales from
126 # http://terraserver.homeadvisor.msn.com/
127 # /About/AboutLinktoHtml.htm
128 #
129 # Fields
130 # Resolution -- Resolution of the
131 # map in meter per pixel
132 # factor -- Scale factor to turn UTM into
133 # tile number
134 # doq -- Aerial photo available
135 # drg -- Topo map available
136 #
137 my %scale_info = (
138 10 => {
139 resolution => 1,
140 factor => 200,
141 doq => 1,
142 drg => 0
143 },
144 11 => {
145 resolution => 2,
146 factor => 400,
147 doq => 1,
148 drg => 1
149 },
150 12 => {
151 resolution => 4,
152 factor => 800,
153 doq => 1,
154 drg => 1
155 },
156 13 => {
157 resolution => 8,
158 factor => 1600,
159 doq => 1,
160 drg => 1
161 },
162 14 => {
163 resolution => 16,
164 factor => 3200,
165 doq => 1,
166 drg => 1
167 },
168 15 => {
169 resolution => 32,
170 factor => 6400,
171 doq => 1,
172 drg => 1
173 },
174 16 => {
175 resolution => 64,
176 factor => 12800,
177 doq => 1,
178 drg => 1
179 },
180 17 => {
181 resolution => 128,
182 factor => 25600,
183 doq => 0,
184 drg => 1
185 },
186 18 => {
187 resolution => 256,
188 factor => 51200,
189 doq => 0,
190 drg => 1
191 },
192 19 => {
193 resolution => 512,
194 factor => 102400,
195 doq => 0,
196 drg => 1
197 }
198 );
199 ################################################
200 # map_to_tiles()
201 #
202 # Turn a map into a set of URLs
203 #
204 # Returns the url array
205 ################################################
206 sub map_to_tiles()
207 {
208 my @result;
209
210 # Get the coordinates as UTM
211 my ($zone,$easting,$north)=latlon_to_utm(
212 'GRS 1980',$center_lat, $center_long);
213
214 # Fix the zone, it must be a number
215 $zone =~ /(d+)/;
216 $zone = $1;
217
218 # Compute the center tile number
219 my $center_x =
220 int($easting /
221 $scale_info{$scale}->{factor});
222
223 my $center_y =
224 int($north /
225 $scale_info{$scale}->{factor});
226
227 # Compute the starting location
228 my $start_x = $center_x - int($x_size / 2);
229 my $start_y = $center_y - int($y_size / 2);
230
231 # Compute the ending location
232 my $end_x = $start_x + $x_size;
233 my $end_y = $start_y + $y_size;
234
235 for (my $y= $end_y-1; $y >= $start_y; --$y) {
236 for (my $x = $start_x;
237 $x < $end_x; ++$x) {
238
239 push (@result, {
240 T => $map_type,
241 S => $scale,
242 X => $x,
243 Y => $y,
244 Z =>$zone}
245 );
246 }
247 }
248 return (@result);
249 }
250
251 ################################################
252 # get_file($) -- Get a photo file from an URL
253 #
254 ################################################
255 sub get_file($)
256 {
257 my $url = shift; # URL to get
258
259 # The name of the file we are going to
260 # write into the cache
261 my $file_spec =
262 "$cache_dir/t=$url->{T}_s=$url->{S}_".
263 "x=$url->{X}_y=$url->{Y}_".
264 "z=$url->{Z}.jpg";
265 if (! -f $file_spec) {
266 # Connection to the remote site
267 my $http = new HTTP::Lite;
268
269 # The image to get
270 my $image_url =
271 "http://terraserver-usa.com/tile.ashx?".
272 "T=$url->{T}&S=$url->{S}&".
273 "X=$url->{X}&Y=$url->{Y}&Z=$url->{Z}";
274 print "Getting $image_url
";
275
276 # The request
277 my $req = $http->request($image_url);
278 if (not defined($req)) {
279 die("Could not get url $image_url");
280 }
281
282 # Dump the data into a file
283 my $data = $http->body();
284 open (OUT_FILE, ">$file_spec") or
285 die("Could not create $file_spec");
286 print OUT_FILE $data;
287 close OUT_FILE;
288 }
289 return ($file_spec);
290 }
291
292 ################################################
293 # toggle_type -- Change the map type
294 ################################################
295 sub toggle_type()
296 {
297 if ($map_type == MAP_TOPO) {
298 if ($scale_info{$scale}->{doq}) {
299 $map_type = MAP_PHOTO;
300 }
301 } else {
302 if ($scale_info{$scale}->{drg}) {
303 $map_type = MAP_TOPO;
304 }
305 }
306 }
307
308 ################################################
309 # get_scale_factor -- Get the current scale
310 ################################################
311 sub get_scale_factor()
312 {
313 return ($scale_info{$scale}->{factor});
314 }
315
316 ################################################
317 # set_map_scale($scale) -- Set the scale for map
318 #
319 # Returns
320 # true if the scale was set,
321 # false if it's not possible to set
322 # the scale to the give value
323 ################################################
324 sub set_map_scale($)
325 {
326 # The scale we want to have
327 my $new_scale = shift;
328
329 if (not defined($scale_info{$new_scale})) {
330 return(0);
331 }
332 if ($map_type == MAP_TOPO) {
333 if (not $scale_info{$new_scale}->{drg}) {
334 return(0);
335 }
336 } else {
337 if (not $scale_info{$new_scale}->{doq}) {
338 return(0);
339 }
340 }
341 $scale = $new_scale;
342 return (1);
343 }
344
345 ################################################
346 # scale_exists($scale)
347 #
348 # Return true if the scale exists for
349 # this type of map
350 #################################################
351 sub scale_exists($)
352 {
353 my $test_scale = shift; # Scale to check
354
355 if ($map_type == MAP_TOPO) {
356 if(not $scale_info{$test_scale}->{drg}) {
357 return (0);
358 }
359 } else {
360 if(not $scale_info{$test_scale}->{doq}) {
361 return (0);
362 }
363 }
364 return (1);
365 }
366 ################################################
367 # get_scales -- Get an array of possible scales
368 ################################################
369 sub get_scales()
370 {
371 return ( sort {$a <=> $b} keys %scale_info);
372 }
373 ################################################
374 # move_map($x, $y) -- Move the map in
375 # the X and Y direction
376 ################################################
377 sub move_map($$)
378 {
379 my $x = shift; # Amount to move in X tiles
380 my $y = shift; # Amount to move in Y tiles
381
382 my ($zone,$east,$north)=
383 latlon_to_utm('GRS 1980',
384 $center_lat, $center_long);
385
386 $east -= $x * get_scale_factor();
387 $north -= $y * get_scale_factor();
388
389 ($center_lat, $center_long) =
390 utm_to_latlon('GRS 1980',
391 $zone, $east, $north);
392 }
393 ################################################
394 # cache_dir -- Return the cache directory
395 ################################################
396 sub cache_dir()
397 {
398 return($cache_dir);
399 }
400 ################################################
401 # init_map -- Init the mapping system.
402 ################################################
403 sub init_map()
404 {
405 if (! -d $cache_dir) {
406 if (not mkdir($cache_dir, 0755)) {
407 die("Could not create cache directory");
408 }
409 }
410 }
411
412 1;
413
Using the Module
The first thing you do is call init_map to initialize the module. The mapping system assumes that you have a 3×3-tile topographical map centered around the Grand Canyon.
At this point, you can call map_to_tiles and get a set of image specifications for this map (nine tiles for your 3×3 map). To turn a specification into a file, call get_file.
The function move_map will move the map a certain number of tiles in any direction. If you want to go to a different place entirely, call set_center_lat_long.
You use the toggle_type function to change from a topographical map to an aerial photograph.
Finally, the scale of the map can be adjusted using set_map_scale.
These are the major pieces; we'll get into some of the nasty details in the section "How It Works."
The USGS is responsible for mapping the nation. The folks there are the ones who produce topographical maps. Microsoft maintains a web server that allows you to download a topographical map or aerial photograph for any place in the United States.
The Results
The result is that when you call map_to_tiles, you pass to get_file to get a set of files that you can put together to make a map.
How It Works
The USGS data is online and can be accessed by anyone. Instructions on how to download this data can be found at:
http://terraserver-usa.com/about.aspx?n=AboutLinktoHtml
Coordinate Systems
Earth is not flat. This is a big problem for mapmakers because maps are flat. Most people locate a point on Earth using longitude and latitude. However, these units suffer from some limitations. For example, it's difficult to measure the distance between two longitudes.
Mapmakers would much rather deal with a flat Earth than a round one. For small patches, it's OK to pretend that Earth is flat. So the standard makers have devised a rectangular coordinate system for mapping points on Earth called the Universal Transverse Mercator (UTM) system. There are several different versions of this coordinate system out there and each one uses its own ellipsoid for coordinate conversion. The United States Geological Survey uses the North American Datum of 1983 (NAD83) version.
Perl has a module to convert longitude/latitude to UTM. But there's a problem. This module has no provision for the NAD83 ellipsoid. Turns out that that NAD83 is the same as an earlier standard, the Geodetic Reference System 1980 (GRS 1980). (It took me about three weeks of searching the Web to discover that GRS 1980 and NAD83 are the same. Obviously, Perl programmers aren't the only ones who can be a bit cryptic.)
Figuring out the language the various mapping agencies are using and all the abbreviations is half the battle. The other half is Perl code.
Downloading Map Tiles
From the TerraServer you can download a 200×200-pixel tile containing a map or aerial photograph of any place in the United States. But you need to know the name of the tile. The first step in the process is to turn the longitude/latitude coordinate into the UTM coordinate used by the server:
210 # Get the coordinates as UTM 211 my ($zone,$easting,$north)=latlon_to_utm( 212 'GRS 1980',$center_lat, $center_long);
To download a tile, you need to know five numbers:
| X | The easting number divided by a scale factor |
| Y | The northing number divided by a scale factor |
| Z | The zone number |
| S | The scale factor |
| T | The map type (1=Topographical, 2= Aerial Photograph, 3=Urban Aerial Photographs) |
Table 10-1 shows the various scale factors for each zoom level.
| Theme | Scale Value | Resolution (Meters per Pixel) | UTM Multiplier |
|---|---|---|---|
| Urban | 8 | 0.25 | 50 |
| Urban | 9 | 0.5 | 100 |
| DOQ, Urban | 10 | 1 | 200 |
| DOQ, DRG, Urban | 11 | 2 | 400 |
| DOQ, DRG, Urban | 12 | 4 | 800 |
| DOQ, DRG, Urban | 13 | 8 | 1,600 |
| DOQ, DRG, Urban | 14 | 16 | 3,200 |
| DOQ, DRG, Urban | 15 | 32 | 6,400 |
| DOQ, DRG, Urban | 16 | 64 | 12,800 |
| DOQ, DRG, Urban | 17 | 128 | 25,600 |
| DOQ, DRG, Urban | 18 | 256 | 51,200 |
| [1] | |||
[1] From the API specification: http://terraserver-usa.com/about.aspx?n=AboutLinktoHtml
The TerraServer contains three types of data. The first, digital raster graphic (DRG), is a topographical map. The next, digital orthophoto quadrangle (DOQ), is an aerial photograph. Finally there is Urban, which indicates a USGS Urban Area photograph. This script does not handle Urban images because they cover only a limited area and because at the time the script was originally written, this type of data was not available.
So let's see what it takes to create a map of the Grand Canyon. You start with the coordinates of the visitor's center in the park:
- 36°03′20″N 112°08′20″W
Now you need to get the S, T, X, Y, and Z values for the tile. You want a topographical map, so the type is 1 (T=1), and you want the highest resolution possible. For topographical maps, that is 1 meter per pixel. Looking through the table, you can see that the scale factor is 11 (S=11).
When you convert the longitude/latitude to UTM, you get this:
| Zone | 12S |
| Easting | 397424 |
| Northing | 3990710 |
The TerraServer wants the zone without the north/south indicator, so the zone is 12 (Z=12).
The table shows that the scale factor is 800. Dividing that into the easting, you get 496 (X=496). Performing a similar conversion on the northing gives you a Y of 4988. As a result, the full URL for the map tile is http://terraserver-usa.com/tile.ashx?T=2&S=12&X=496&Y=4988&Z=12.
Note:
The X- and Y-coordinate numbers are consecutive. So by decrementing the X number by 1, you get the tile to the left of the current tile, incrementing the Y number by 1 gives the tile below the current tile, and so on.
Getting the Data
The get_file function is responsible for turning a tile specification into an image file. The module HTTP::Lite is used to fetch the file.
The first thing you do is create a HTTP::Lite object for downloading:
266 # Connection to the remote site 267 my $http = new HTTP::Lite;
Next you turn your tile specification into a URL:
269 # The image to get
270 my $image_url =
271 "http://terraserver-usa.com/tile.ashx?".
272 "T=$url->{T}&S=$url->{S}&".
273 "X=$url->{X}&Y=$url->{Y}&Z=$url->{Z}";
The next step is to create an HTTP request to get the data:
276 # The request
277 my $req = $http->request($image_url);
278 if (not defined($req)) {
279 die("Could not get url $image_url");
280 }
This gets all sorts of information about the page. All you want is the data, so you take the body of the page and dump it to a file. It is this file that you give back to the user as the image file they want:
282 # Dump the data into a file
283 my $data = $http->body();
284 open (OUT_FILE, ">$file_spec") or
285 die("Could not create $file_spec");
286 print OUT_FILE $data;
287 close OUT_FILE;
288 }
Moving the Map
You allow the map to be panned to the left or right. The move_map function moves the map by tiles. But you store your center point as longitude/latitude. Changing the center is not as simple as just adding in a constant to these values.
The problem is that longitude curves. So in order to recenter, you need a rectangular coordinate system, in this case UTM. The amount to move is determined by the scale factor. The move_map function schanges the center point by one tile in the X or Y direction or both. Each parameter to this function can have the value 1, 0, or -1. The result of this function is a new map with a different center point.
373 ################################################
374 # move_map($x, $y) -- Move the map in
375 # the X and Y direction
376 ################################################
377 sub move_map($$)
378 {
379 my $x = shift; # Amount to move in X tiles
380 my $y = shift; # Amount to move in Y tiles
381
382 my ($zone,$east,$north)=
383 latlon_to_utm('GRS 1980',
384 $center_lat, $center_long);
385
386 $east -= $x * get_scale_factor();
387 $north -= $y * get_scale_factor();
388
389 ($center_lat, $center_long) =
390 utm_to_latlon('GRS 1980',
391 $zone, $east, $north);
392 }
Hacking the Script
This module was created by the process of successive experimentation: try something, see if works, try something else, see if it works, add a little to the code, and so on. In other words, there's not a whole lot of design that went into this module.
As a result, the API is a little more complex and cluttered than it needs to be. The code could use a little cleaning up. But then again, this is Wicked Cool Perl Scripts, not Clean Pretty Perl Scripts, so have fun.