#!/usr/bin/perl -w ##################################################################### # # Program name: genData.pl # # Description : script to fetch earthquake bulletin # and convert to shapefile # # $Id: genData.pl,v 1.6 2004/01/02 21:20:30 tkralidi Exp $ # # Revisions: (see end of file for revision history) # ##################################################################### # enable strict mode use strict; # load URL fetching module use LWP::Simple; # load shapefile creation module use mapscript; # load DBase creation module use XBase; # set variables my @quakes; my $i = 0; my $basedir = "/home/tkralidi/neis/"; my $basename = $basedir . "neis"; # URL to raw data my $url = "http://neic.usgs.gov/neis/finger/quake.asc"; # fetch data and store in array my $content = get($url); my @arr = split /\n/, $content; # loop through file and pick up attributes # we have to put a few conditions here # to realize empty data foreach my $line(@arr) { my ($date, $time, $x, $y); if ($line =~ /^0/) { # pick up attribute info # if "Q" is empty if ($line =~ /^(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s{5}(.*)$/) { $date = $1; $time = $2; $y = $3; $x = $4; $quakes[$i]{depth_km} = $5; $quakes[$i]{mag} = $6; $quakes[$i]{q} = "none"; $quakes[$i]{comment} = $7; } # if "MAG" and "Q" are empty elsif ($line =~ /^(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s{10}(.*)$/) { $date = $1; $time = $2; $y = $3; $x = $4; $quakes[$i]{depth_km} = $5; $quakes[$i]{mag} = 9999; # dummy numerical value $quakes[$i]{q} = "none"; $quakes[$i]{comment} = $6; } # elsif "MAG" and "Q" are NOT empty elsif ($line =~ /^(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(.*)$/) { $date = $1; $time = $2; $y = $3; $x = $4; $quakes[$i]{depth_km} = $5; $quakes[$i]{mag} = $6; $quakes[$i]{q} = $7; $quakes[$i]{comment} = $8; } # alter coords to be expressed as +/- vals if ($y =~ /^(\d+\.\d+)N$/) { $y = $1; } if ($y =~ /^(\d+\.\d+)S$/) { $y = "-" . $1; } if ($x =~ /^(\d+\.\d+)E$/) { $x = $1; } if ($x =~ /^(\d+\.\d+)W$/) { $x = "-" . $1; } $quakes[$i]{y} = $y; $quakes[$i]{x} = $x; # format ISO8160 timestring like YYYY-MM-DDTHH:MM:SSZ $date =~ s/\//-/g; $quakes[$i]{datetime} = "20" . $date . "T" . $time . "Z"; # get rid of magnitude unit symbol 'M' so we can classify # with numerical comparisons $quakes[$i]{mag} =~ s/M//; $i++; } } # now that we have the data, remove previous copies of data unlink("$basename.shp"); unlink("$basename.shx"); unlink("$basename.dbf"); # make new .shp/.shx/.dbf my $newshpf = new mapscript::shapefileObj($basename, $mapscript::MS_SHAPEFILE_POINT); my $newtable = XBase->create("name" => $basename, "field_names" => [ "DATETIME", "DEPTH_KM", "MAG", "Q", "COMMENT" ], "field_types" => [ "C", "N", "N", "C", "C" ], "field_lengths" => [ 20, 10, 10, 5, 100 ], "field_decimals" => [ undef, 5, 5, undef, undef ]); $i = 0; # populate records foreach my $k (@quakes) { # add to .dbf my @data = $newtable->set_record($i, $k->{datetime}, $k->{depth_km}, $k->{mag}, $k->{q}, $k->{comment}); # add to .shp/.shx my $point = new mapscript::pointObj(); $point->{x} = $k->{x}; $point->{y} = $k->{y}; $newshpf->addPoint($point); print "$k->{datetime} -- $k->{y} -- $k->{x} -- $k->{depth_km} -- $k->{mag} -- $k->{q} -- $k->{comment}\n"; $i++; } ##################################################################### # # Revision History # # $Log: genData.pl,v $ # #####################################################################