Skip to content

anonymous /yearofedits.py

Embed URL

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Given a OpenStreetMap changeset download, and a population density download, plots the global edits-per-population density
#!/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)
@danstowell

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Something went wrong with that request. Please try again.