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]:
[<matplotlib.lines.Line2D at 0x110834250>]