mapquakes.pl

Here's the pilot script that MapScript was written from. It's called mapquakes.pl and is an example that integrates remote data (earthquake locations from the USGS) with a local raster image of the earth. Here's how it works:

1) 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

2) parse the file, storing the location and magnitude information

3) draw the raster layer (global topographic data)

4) 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.

5) 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

---

#!/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