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.