Turning GSHAP data into a geoImage

README.md

Data GSHAP data to image

This is the data downloading and cleaing parts of the map making process described in this blog post.

to run, open your terminal

$ make

create_earthquake_raster.py

#!/usr/bin/env python
# encoding: utf-8

import Image

def main():
	#open the data file
	with open("GSHPUB.DAT") as f:
		data = f.read().split("\r\n")
	data = map(lambda x: x.split("\t"), data)

	#define the dimenstions of the image
	#one degree for each 1/10th of a degree
	width = 180*2*10
	height = 90*2*10

	#create a new image object
	img = Image.new("RGB",[width+1,height+1])
	pixels = img.load()

	#loop through the data
	for d in data:
		try:
			#shift the data to be all positive numbers
			x = int(float(d[0])*10) + (180*10)
			y = height - (int(float(d[1])*10) + (90*10))

			#set no data to 0
			c = int(float(d[2])*10) if d[2] != "NaN" else 0

			#if the data is set too far east...wrap it to the west
			if x > width:
				x = x - width

			#assign the data and check for errors
			try:
				pixels[x,y] = (c,c,c)
			except IndexError:
				print "Error with: ", x,y,c
		except ValueError:
			print "Error with: ", d

	#save the image
	img.save("output.tiff", "TIFF")
	
	pass

if __name__ == '__main__':
	main()

Makefile

DATA_FILE = ./GSHPUB.DAT
DATA_ZIP = ./gshpub.zip
DATA_ZIP_LOCATION = http://www.seismo.ethz.ch/static/gshap/gshpub.zip

IMAGE_SCRIPT = ./create_earthquake_raster.py

OUTPUT_IMAGE = ./output.tiff

all: ${OUTPUT_IMAGE}


${OUTPUT_IMAGE}: ${IMAGE_SCRIPT} ${DATA_FILE}
	@echo creating image file
	python ${IMAGE_SCRIPT}

${DATA_FILE}: ${DATA_ZIP}
	@echo unzipping data file
	unzip -o ${DATA_ZIP}

${DATA_ZIP}:
	@echo downloading data file
	curl ${DATA_ZIP_LOCATION} -O -L


clean:
	rm ${DATA_FILE}
	rm ${OUTPUT_IMAGE}
	rm ${DATA_ZIP}