In [37]:
import json
from math import pi, log, tan, ceil
Function to convert lng, lat from geojson to mercator XY
In [108]:
def lnglatToXY(ll):
D2R = pi / 180.0
A = 6378137.0
x = (A * ll[0] * D2R)
y = (A * log(tan((pi*0.25) + (0.5 * ll[1] * D2R))))
return x,y
Read GeoJSON, insert projected coords into array
In [23]:
with open('/Users/dnomadb/old-make-surface/python-processors/sf_bike_crimes.geojson', 'r') as pointFile:
geojson = json.loads(pointFile.read())
xys = []
for f in geojson['features']:
xys.append(lnglatToXY(f['geometry']['coordinates']))
del geojson
xys = np.array(xys)
Cellsize = dimensions, in meters, of the analysis grid
In [109]:
cellsize = 1000
In [12]:
xmin = xys[:,0].min()
xmax = xys[:,0].max()
ymin = xys[:,1].min()
ymax = xys[:,1].max()
x_range = xmax - xmin
y_range = ymax - ymin
Estimate the row number, col number, and a few other things for processing
In [110]:
cols = int(ceil(x_range / float(cellsize)))
rows = int(ceil(y_range / float(cellsize)))
colInds = arange(cols + 1)
rowInds = arange(rows + 1)
colBins = colInds * cellsize + x_origin
rowBins = rowInds * cellsize + y_origin
In [119]:
fig = figure(figsize=(11,10))
rnded = np.round(xys / float(cellsize)) * cellsize
hist = np.histogram2d(xys[:,0], xys[:,1], [colBins, rowBins])
for r in range(rows):
for c in range(cols):
shparr = array([
[hist[1][c], hist[2][r]],
[hist[1][c], hist[2][r + 1]],
[hist[1][c + 1], hist[2][r + 1]],
[hist[1][c + 1], hist[2][r]],
[hist[1][c], hist[2][r]]
])
## array of indices that are in this particular bin
geoJSONindices = where((rnded[:,0] == hist[1][c]) & (rnded[:,1] == hist[2][r]))[0]
## eg to see number of crimes in the cell:
## print len(geoJSONindices),
## plot to see the grid
plot(shparr[:,0], shparr[:,1], c='grey')
## plot the points to see where they are
plot(xys[:,0], xys[:,1], '+')
Out[119]: