createkml.py
A python script that creates a kml representation of earthquake data.
createkml.py expects a path to a csv file of the form created by the Geonet csv export call.
The implementation I use to do this reads the file downloaded by a bash script (running as a cron job)
The essential part of this is
wget -O /PATH/TO/FILE.csv "http://wfs.geonet.org.nz/geonet/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=geonet:quake_search_v1&maxFeatures=50&outputFormat=csv"
The format of the file has the following columns, we are only interested in a few of them (in bold)
FID,
publicid,
eventtype,
origintime,
modificationtime,
latitude,
longitude,
depth,
magnitude,
evaluationmethod,
evaluationstatus,
evaluationmode,
earthmodel,
depthtype,
originerror,
usedphasecount,
usedstationcount,
minimumdistance,
azimuthalgap,
magnitudetype,
magnitudeuncertainty,
magnitudestationcount,
origin_geom
See http://info.geonet.org.nz/
The code
Basic stuff, no data integrity checking (this is part of a pipeline in the local implementation).
#!/usr/bin/python import argparse import StringIO import csv import datetime import dateutil.parser import math # build the parser parser = argparse.ArgumentParser(description='Plot earthquakes on a map') parser.add_argument('filename', help='path/to/csv file') args = parser.parse_args() # open the file provided f = open(args.filename, 'r') # read it in, and create the basic data we will # be using output = StringIO.StringIO(f.read()) # import in csv file format data = csv.reader(output) # drop the first line data.next() # print out the top of the file print '<!--?xml version="1.0" encoding="UTF-8"?-->\n\n' print '\tRecent earthquakes in and around New Zealand' print '\tSource: GeoNet (www.geonet.org.nz); Created: {}'.format(repr(datetime.datetime.now().isoformat())) # get the time right noe (UTC) now = datetime.datetime.utcnow() # maxdepth is the deepest value # after which all colour values are clamped maxdepth = 400.0 for row in data: # convert the date to something humans understand whendatetime = dateutil.parser.parse(row[3]) # figure out how far to fade the icon # clamp this to ~25% minimum difference = (now - whendatetime).total_seconds(); if difference > 86400: difference = 86400 elif difference < 0: difference = 0 # 'shade' is a two character hex value, 3f-ff, # used as the alpha value in the colour applied # to the icon shade = format(255 - int((float(difference)/86400)* 192),'02x') # Generate a colour value for the icon # associated with depth. # The color ramp will be a logarithmic # progression from 00ff00 (shallow) # to 0000ff (deep), maxed at maxdepth km? depth = float(row[7]) # fast at shallow, slow at depth # multiplier = 255/sqrt(maxdepth) # depth2 = math.sqrt(depth)*(255/math.sqrt(maxdepth)) # slow at shallow, fast at depth # divisor = maxdepth/sqrt(255) # depth2 = (depth/(maxdepth/math.sqrt(255)))**2 # flat # multiplier = maxdepth/255 depth2 = depth*(255/maxdepth) if depth2 > 255: depth2 = 255.0 elif depth2 < 0: depth2 = 0.0 if depth2 < 128: green = 'ff' red = format(int(depth2*2),'02x') else: green = format(255-int((depth2-128)*2),'02x') red = 'ff' # round magnitude to determine the icon style magnitude = float(row[8]); magnitudelow = int(math.floor(magnitude)) if magnitudelow > 9: magnitudelow = 9 elif magnitudelow < 0: magnitudelow = 0 description = '<![CDATA[\n <div style="text-align: left;">' description += '<a href="http://www.geonet.org.nz/quakes/region/newzealand/' + row[1] + '" target="' + row[1] + '"><b>' + row[1] + '</b></a>' description += '\n <b>' + whendatetime.strftime('%H:%M, %d %B %Y') + '</b> (UTC)' description += '\n magnitude: <b>' + str(round(float(row[8]),2)) + '</b>' description += '\n depth: <b>' + str(round(float(row[7]),2)) + 'km</b>' description += '</div> ]]>' print '\t' print '\t\t' + description + '' print '\t\t' print '\t\t\n\t\t\t' + row[6] + ',' + row[5] + ',0' print '\t\t\n\t' # print out the bottom of the file print '\n\n'
The script dumps its output to standard output, so I generally redirect this to a file.