Given a OpenStreetMap changeset download, and a population density download, plots the global edits-per-population density

yearofedits.py

#!/usr/bin/env python

from xml.sax import make_parser
from xml.sax.handler import ContentHandler
from math import floor, log
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from colorsys import hsv_to_rgb
import numpy as np
import os.path
import csv

##################
whichyear = 2012
gridw = 360 * 2
gridh = 180 * 2

##################
# load the GPW population-density data
popgridfile = 'gpw/glds00ag30.asc'
rdr = csv.reader(file(popgridfile, 'rb'), delimiter=' ')
# handle header lines
pophdr = {}
for _ in range(6):
	row = rdr.next()
	if row[0] in ['NODATA_value', 'cellsize']:
		pophdr[row[0]] = float(row[-1])
	else:
		pophdr[row[0]] =   int(row[-1])

# now the data
popgrid = [[float(datum) for datum in row[:-1]] for row in rdr]

if np.shape(popgrid) != (pophdr['nrows'], pophdr['ncols']):
	raise ValueError("popdata shape mismatch against header declaration: %s != %s" % (np.shape(popdata), (pophdr['nrows'], pophdr['ncols'])))
print "popgrid shape:"
print np.shape(popgrid)

popgrid_log = [[log(max(1e-6, datum)) for datum in row] for row in popgrid]

plt.figure()
plt.imshow((popgrid_log),
	    cmap=cm.gist_earth_r,
	    aspect='auto')
plt.title("Population density", fontsize=10)
plt.xticks([])
plt.yticks([])
plt.savefig('plot_populationdensity.png', papertype='A4', format='png')

##################
class ChangesetsHandler(ContentHandler):
	def __init__ (self, year):
		self.year = year
		self.yearstr = str(year)
		self.gridw = gridw
		self.gridh = gridh
		self.grid = [[0.0 for _ in range(gridh)] for _ in range(gridw)]
		print "grid shape:"
		print np.shape(self.grid)
		self.lonadd = 180
		self.lonmul = float(gridw) / 360.0
		self.latadd = 90
		self.latmul = float(gridh) / 180.0
	def startElement(self, name, attrs):
		if (name == 'changeset') and (attrs.get('num_changes', '0')!='0'):
			if (attrs.get('created_at', "")[:4] == self.yearstr):
				ch = {}
				for getit in ['id', 'created_at', 'min_lon', 'min_lat', 'max_lon', 'max_lat']:
					ch[getit] = attrs.get(getit, "")
					if ch[getit]=='':
						#print "warning: %s empty. id: %s, num_changes: %s" % (getit, attrs.get('id'), attrs.get('num_changes'))
						return
				# Here we calculate its bbox in relation to the grid
				x1 = int(floor((float(ch['min_lon']) + self.lonadd) * self.lonmul))
				x2 = int(floor((float(ch['max_lon']) + self.lonadd) * self.lonmul))
				y1 = int(floor((float(ch['min_lat']) + self.latadd) * self.latmul))
				y2 = int(floor((float(ch['max_lat']) + self.latadd) * self.latmul))
				x2 = max(x2, x1+1)
				y2 = max(y2, y1+1)
				density = 1.0 / ((x2-x1)*(y2-y1))
				for x in range(x1, x2):
					for y in range(y1, y2):
						try:
							self.grid[x][y] += density
						except:
							pass

def ratiovalcalc(editval, popval):
	"Calculates a pixel value for the ratiogrid"
	if (popval <= 0.0) or (editval==0.0):
		return -15  #float('nan') # 0.0
	else:
		return log(editval / popval)

def compare_edits_to_pop(editsgrid, popgrid, pophdr):
	editsgrid = np.rot90(np.array(editsgrid))
	popgrid = np.array(popgrid)
	ystart = 10  # MAGIC NUMBER

	width  = np.shape(popgrid)[1]
	height = np.shape(popgrid)[0]

	# shrink editsgrid to match the popgrid
	editsgrid = editsgrid[ystart:ystart+height,:]

	# Manually find the min and max (since the sea makes it tricky...)
	pixelmin = 9999999999
	pixelmax = -9999999999
	pixelavg = 0.0
	for y in range(height):
		for x in range(width):
			if (popgrid[y,x] > 0.0) and (editsgrid[y,x]!=0.0):
				datum = ratiovalcalc(editsgrid[y,x], popgrid[y,x])
				pixelavg += datum
				if pixelmin > datum:
					pixelmin = datum
				if pixelmax < datum:
					pixelmax = datum
	pixelavg /= (height * width)

	# What we want is for the 0--1 range to be the bigger of the above-and-below ranges, so:
	pixelscaler = 1.0 / max(pixelmax-pixelavg, pixelavg-pixelmin)
	# but we also exaggerate a bit to emphasise the middle of the range:
	pixelscaler *= 2

	# Now we create the RGB values
	ratiogrid = [[0 for _ in range(width)] for _ in range(height)]
	for y in range(height):
		for x in range(width):
			if (popgrid[y,x] > 0.0) and (editsgrid[y,x]!=0.0):
				datum = ratiovalcalc(editsgrid[y,x], popgrid[y,x])
				#datum = ((datum - pixelmin) * 2 / (pixelmax - pixelmin)) - 1
				datum = (datum - pixelavg) * pixelscaler
				if datum > 0:
					colour = hsv_to_rgb(0.6, max(min( datum, 1.0), -1.0), 1)
				else:
					colour = hsv_to_rgb(0.0, max(min(-datum, 1.0), -1.0), 1)
			else:
				colour = [0,0,0]
			ratiogrid[y][x] = colour

	plt.imshow(ratiogrid,
		    aspect='auto')
	plt.title("Relative edits per unit population (log scale), %i" % whichyear, fontsize=10)
	plt.xticks([])
	plt.yticks([])
	plt.savefig('plot_edits_per_pop.png', papertype='A4', format='png')

if __name__ == '__main__':
	parser = make_parser()   
	handler = ChangesetsHandler(whichyear)
	parser.setContentHandler(handler)
	parser.parse(open('/home/dan/dev/osm/changesets/changesets-latest.2012only.osm'))

	print "Shape of grid got:"
	print np.shape(handler.grid)

	# loggify
	editsgrid = handler.grid
	editsgrid_log = [[log(max(1e-6, datum)) for datum in row] for row in editsgrid]

	plt.figure()
	plt.imshow(np.rot90(editsgrid_log),
		    cmap=cm.gist_earth_r,
		    aspect='auto')
	plt.title("Distribution of OpenStreetMap edits (bounding boxes), %i" % whichyear, fontsize=10)
	plt.xticks([])
	plt.yticks([])

	plt.savefig('plot_yearofedits_%i.png' % whichyear, papertype='A4', format='png')

	compare_edits_to_pop(editsgrid, popgrid, pophdr)

this is by me. seems I wasn't logged in when I pasted.