public
anonymous / yearofedits.py
Created

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

  • Download Gist
yearofedits.py
Python
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172
#!/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.

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.