forked from MapServer/MapServer
-
Notifications
You must be signed in to change notification settings - Fork 2
PerlMapScriptExamples35ex2
Thomas Bonfort edited this page Apr 6, 2012
·
2 revisions
-
use the perl LWP module to retrieve a text file from USGS servers that contains the locations and magnitudes of the most recent earthquakes worldwide
-
parse the file, storing the location and magnitude information
-
draw the raster layer (global topographic data)
-
loop through the quakes drawing each in turn, scaling the symbol by the magnitude. Each quake is labeled with an ID that is used to reference the quake in a lookup table that is presented below the interface.
-
save the final map, and return a page to the user
The script shows how to build a simple pan/zoom interface and details how one might build maps from disparate sources. Map files and topographic data are available upon request. This should run just fine with !MapServer version 3.5 and I will post a version that runs with 3.6 soon.
-- Steve
#!perl
#!/usr/bin/perl
use CGI qw(:standard :html escape);
use LWP::Simple;
use mapscript;
#
# defaults, make any necessary changes here...
#
$quake_map = '/home/sdlime/mapquakes/quakes.map';
$quake_source = 'ftp://ghtftp.cr.usgs.gov/pub/quake';
@quakes = '';
$quake_data = '';
$quake_timestamp = '';
$image_path = '/usr/local/apache/htdocs/tmp/';
$image_virtual_path = '/tmp/';
$image_id = $$ . time() . ".gif";
$map = '';
$zoom_direction = 0;
$zoom_size = 2;
sub get_quakes {
$quake_data = get $quake_source or die 'get failed';
my $i = 0;
my $is_quake = 0;
foreach (split(/^/, $quake_data)) {
$quake_timestamp = $_ if /^Updated as of/;
if($is_quake) {
$quakes[$i]{lat} = substr($_, 18, 7);
$quakes[$i]{lon} = substr($_, 26, 7);
$quakes[$i]{magnitude} = substr($_, 40, 3);
if($quakes[$i]{lon} =~ /E$/) {
$quakes[$i]{lon} =~ s/E//;
} else {
$quakes[$i]{lon} =~ s/W//;
$quakes[$i]{lon} *= -1;
}
if($quakes[$i]{lat} =~ /N$/) {
$quakes[$i]{lat} =~ s/N//;
} else {
$quakes[$i]{lat} =~ s/S//;
$quakes[$i]{lat} *= -1;
}
$i++;
}
$is_quake = 1 if /^yy/;
}
return;
}
sub set_extent() {
my $zoom, @imgext;
my $x, $y;
my $cellsizex, $cellsizey;
if($cgi->param('imgext')) { # if an interactive interface then calculate a new extent
@imgext = split(' ', $cgi->param('imgext'));
$x = $cgi->param('img.x');
$y = $cgi->param('img.y');
$zoom_size = $cgi->param('zoomsize');
$zoom_direction = $cgi->param('zoomdir');
if($zoom_direction == 0) { # pan
$zoom = 1;
} else { # in or out
$zoom = $zoom_size*$zoom_direction;
$zoom = 1.0/abs($zoom) if $zoom < 0;
}
$cx = ($imgext[2]-$imgext[0])/($map->{width}); # calculate cellsize in x and y
$cy = ($imgext[3]-$imgext[1])/($map->{height});
$x = $imgext[0] + $cx*$x; # change x,y from image to map coordinates, offset from UL corner of previous image
$y = $imgext[3] - $cy*$y;
$map->{extent}->{minx} = $x - .5*(($imgext[2] - $imgext[0])/$zoom); # calculate new extent
$map->{extent}->{miny} = $y - .5*(($imgext[3] - $imgext[1])/$zoom);
$map->{extent}->{maxx} = $x + .5*(($imgext[2] - $imgext[0])/$zoom);
$map->{extent}->{maxy} = $y + .5*(($imgext[3] - $imgext[1])/$zoom);
}
}
sub render_quakes {
die $mapscript::ms_error->{message} unless $map = new mapObj($quake_map);
&set_extent();
my $img = $map->prepareImage();
my $point = new pointObj();
my $layer = $map->getLayerByName('relief');
$layer->draw($map, $img); # draw basemap
$layer = $map->getLayerByName('quakes');
my $class = $layer->getClass(0);
my $i = 1; # quake counter
foreach my $q (@quakes) {
$class->{size} = int($q->{magnitude}*2);
$class->{maxsize} = $class->{minsize} = $class->{size}; # don't want to scale
$point->{x} = $q->{lon};
$point->{y} = $q->{lat};
$point->draw($map, $layer, $img, undef, "$i");
$i++;
}
#$layer = $map->getLayerByName('timestamp'); # add timestamp
#$layer->{status} = 1; # turn the layer on
#$class = $layer->getClass(0);
#$class->setText($quake_timestamp);
#$layer->draw($map, $img);
$map->drawLabelCache($img);
mapscript::msSaveImage($img, $image_path . $image_id, $map->{imagetype}, $map->{transparent}, $map->{interlace}, $map->{imagequality});
return;
}
# here is the main program
$cgi = new CGI;
&get_quakes();
&render_quakes();
# now the html
print header();
print start_html(-title=>'MapServer Test Suite - Earthquake Mapper', -bgcolor=>"#ffffff");
print "<!-- " .mapscript::msGetVersion(). " -->";
print "<!-- " .$map->{scale}. " -->";
print "<center><table width=\"600\"><tr><td>";
print "<font size=\"+2\" face=\"arial,helvetica\"><b>MapScript - Earthquake Mapper</b></font><p>";
print "This is a simple of example of using MapScript. Near real-time earthquake data is retrieved ";
print "from the USGS web servers. The resulting information is parsed using perl and each ";
print "quake is plotted on top of a basemap. The size of the marker indicates quake magnitude. ";
print "Many other possible enhancements. If the interface seems a bit slow it's because the ";
print "USGS servers can be difficult to reach- the perils of remote data. Here's the ";
print "<a href=\"/mapserver_demos/tests/perl/quakes/mapquakes.pl\">source</a> code.<p>";
print "<form name=\"mapquakes\" action=\"/cgi-bin/mapquakes.pl\" method=\"get\">";
print "<center>";
print "<input border=\"0\" type=\"image\" name=\"img\" src=\"". $image_virtual_path . $image_id ."\"><br>";
print "<input type=\"hidden\" name=\"imgext\" value=\"" . join(' ', $map->{extent}->{minx},$map->{extent}->{miny},$map->{extent}->{maxx},$map->{extent}->{maxy}) ."\">";
if($zoom_direction == 0) {
print "<input type=\"radio\" name=\"zoomdir\" value=\"1\"> zoom in ";
print "<input type=\"radio\" name=\"zoomdir\" value=\"0\" checked> pan ";
print "<input type=\"radio\" name=\"zoomdir\" value=\"-1\"> zoom out ";
} else {
if($zoom_direction == -1) {
print "<input type=\"radio\" name=\"zoomdir\" value=\"1\"> zoom in ";
print "<input type=\"radio\" name=\"zoomdir\" value=\"0\"> pan ";
print "<input type=\"radio\" name=\"zoomdir\" value=\"-1\" checked> zoom out ";
} else {
print "<input type=\"radio\" name=\"zoomdir\" value=\"1\" checked> zoom in ";
print "<input type=\"radio\" name=\"zoomdir\" value=\"0\"> pan ";
print "<input type=\"radio\" name=\"zoomdir\" value=\"-1\"> zoom out ";
}
}
print "zoom size <input type=\"text\" size=\"2\" name=\"zoomsize\" value=\"$zoom_size\">";
print "</center>";
print "</form>";
print "<p><font face=\"arial,helvetica\"><b>Here's some of the output from the USGS:</b></font><p>";
print "<p><pre>";
$show = 0;
$number = 0;
foreach (split(/^/, $quake_data)) {
$show = 1 if /Updated/;
printf "<b>%2d</b>: ", $number if $show and $number;
print " ". $_ if $show and !$number;
print $_ if $show and $number;
$number++ if $show and $number;
$number = 1 if /^yy/;
}
print "<\/pre>;";
print "</td></tr></table>";
print end_html();
exit(0);
== Begin Quakes.map ==
MAP
STATUS ON
EXTENT -180 -90 180 90
SIZE 561 281
UNITS DD
IMAGECOLOR 255 0 0
LAYER
NAME relief
TYPE RASTER
DATA relief.tif
STATUS DEFAULT
END
SYMBOL
TYPE ELLIPSE
NAME "circle"
POINTS 1 1 END
FILLED
END
LAYER
NAME quakes
TYPE POINT
STATUS OFF
SYMBOLSCALE 50000000
CLASS
NAME "Recent earth quakes."
SYMBOL circle
SIZE 15
OUTLINECOLOR 0 0 0
COLOR 255 0 0
LABEL
OUTLINECOLOR 0 0 0
COLOR 255 255 255
SIZE TINY
POSITION AUTO
PARTIALS FALSE
END
END
TEMPLATE "dummy"
END
LAYER
NAME timestamp
TYPE ANNOTATION
STATUS OFF
TRANSFORM FALSE
CLASS
TEXT 'This is where the timestamp goes.'
LABEL
COLOR 0 0 0
OUTLINECOLOR 225 225 225
SIZE MEDIUM
POSITION CR
FORCE ON
END
END
FEATURE
POINTS 5 270 END
END
END
END
}}