In [1]:
import shapefile as shp
import math
import pandas as pd
import geopandas as gpd
import pylab as pl
from fiona.crs import from_epsg
import branca.colormap as cm
import mplleaflet
import folium
%pylab inline
In [2]:
# grid
delta = 0
minx,maxx,miny,maxy = 982258.382157, 984756.445048, 209697.251114, 212819.601834
dx = 50
dy = 50
nx = int(math.ceil(abs(maxx - minx)/dx))
ny = int(math.ceil(abs(maxy - miny)/dy))
w = shp.Writer(shp.POLYGON)
w.autoBalance = 1
w.field("ID")
id=0
for i in range(ny):
for j in range(nx):
id+=1
vertices = []
parts = []
vertices.append([min(minx+dx*j,maxx),max(maxy-dy*i,miny)])
vertices.append([min(minx+dx*(j+1),maxx),max(maxy-dy*i,miny)])
vertices.append([min(minx+dx*(j+1),maxx),max(maxy-dy*(i+1),miny)])
vertices.append([min(minx+dx*j,maxx),max(maxy-dy*(i+1),miny)])
parts.append(vertices)
w.poly(parts)
w.record(id)
w.save('polygon_grid')
In [3]:
# read files
grid = gpd.read_file('./polygon_grid.shp')
grid.crs = from_epsg(2263)
#grid.to_crs(epsg=2263, inplace=True)
uni_bssid = pd.read_pickle("./unique_bssid.p")
uni_bssid = gpd.GeoDataFrame(uni_bssid)
uni_bssid.crs = from_epsg(2263)
# uni_bssid.to_crs(epsg=2263, inplace=True)
PointInPoly = gpd.sjoin(uni_bssid, grid, how='left', op='intersects')
# Somehow, three points do not intersect with any grid cell.
# print PointInPoly.count()
# Delete ID= 'NA'
PointInPoly.dropna(subset=['ID'], inplace=True)
# print PointInPoly.count()
PointInPoly.head(2)
# groupby
grouped = PointInPoly.groupby('ID').apply(lambda x: reduce(lambda x,y: x+y, x.unique_bssid_list))
bssidInPoly = pd.DataFrame(grouped, columns=['all_bssid_list'])
bssidInPoly['unique_bssid_list'] = map(lambda x: set(x), grouped)
bssidInPoly['cum'] = map(lambda x: len(x), grouped)
bssidInPoly['uni'] = map(lambda x: len(set(x)), grouped)
bssidInPoly['ID'] = bssidInPoly.index
bssidInPoly.reset_index(drop=True, inplace=True)
#bssidInPoly.head(2)
# merge grid and bssidInPoly
grid_bssid = pd.merge(grid, bssidInPoly, how='left', on='ID')
# grid_bssid.uni.fillna(inplace=True, value=0)
grid_bssid.to_crs(epsg=2263, inplace=True)
grid_bssid.sort_values(by='uni', ascending=False).head(3)
Out[3]:
In [7]:
# Data for Plot
grid_plot = grid_bssid.loc[:, ['ID', 'uni', 'geometry']]
grid_plot.dropna(subset=['uni'], inplace=True)
#grid_plot.uni.fillna(inplace=True, value=0)
# Select colormap
# import branca.colormap as cm
linear = cm.LinearColormap(['#000000', 'red'])
linear
# Colormap setting
colormap = linear.to_step(8)
colormap.caption = '# of bssid (unique)'
# Add 'style' - scale 'uni' values
scale = grid_plot.uni.max() - grid_plot.uni.min()
uni_scale = grid_plot.uni / scale
grid_plot['style'] = [{'fillColor': colormap(i), 'weight': 0, 'color': 'black', 'dashArray': '0, 0'} for i in uni_scale]
In [8]:
m = folium.Map([40.743, -74], zoom_start=15, tiles='Cartodb dark_matter', crs='EPSG3857')
folium.GeoJson(grid_plot, overlay=True).add_to(m)
m.add_child(colormap)
m
Out[8]:
In [9]:
kw = dict(column='uni', k=30, colormap='OrRd')
f, (ax1, ax2, ax3) = pl.subplots(1,3,figsize=(20,5))
d_ax = {'ax1':ax1, 'ax2':ax2, 'ax3':ax3}
d_sch = {'ax1':'QUANTILES', 'ax2':'equal_interval', 'ax3':'fisher_jenks'}
for i in range(1,4,1):
axi = 'ax'+str(i)
grid_plot.plot(scheme=d_sch[axi], ax=d_ax[axi], **kw)
d_ax[axi].set_title(d_sch[axi])
In [11]:
# fig, ax = pl.subplots(1,1,figsize==(5,5))
# grid_plot.plot(ax=ax, colormap='CMRmap', c=grid_plot.uni)
# mplleaflet.display(fig=ax.figure, crs=df.crs)
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [12]:
hp_target = pd.read_pickle("./hp_target.p")
hp_target.crs = from_epsg(2263)
hp_target['style'] = [{'fillColor': '#000000', 'weight': 1, 'color': 'tan'}] * len(hp_target)
In [17]:
hp_target.drop(['style', 'tds_num'], axis=1).to_file('hp_target')
grid_plot.drop(['style'],axis=1).to_file("grid_plot")
In [22]:
grid_bssid.uni.fillna(inplace=True, value=0)
grid_bssid.drop([u'all_bssid_list', u'unique_bssid_list', u'cum'], axis=1).to_file("grid_bssid")
In [13]:
m = folium.Map([40.743, -74], zoom_start=15, tiles='Cartodb dark_matter', crs='EPSG3857')
folium.GeoJson(grid_plot, overlay=False).add_to(m)
folium.GeoJson(hp_target, overlay=True).add_to(m)
m.add_child(colormap)
m
Out[13]:
In [23]:
from IPython.display import IFrame
IFrame('https://nyu.carto.com/u/alanf/builder/29ad8d2e-579b-11e7-af33-0e3a376473ab/embed',
width=900, height=750)
Out[23]:
In [ ]:
In [ ]:
In [53]:
free = gpd.read_file("./hp_plot/hp_plot.shp")
free.crs = from_epsg(2263)
In [54]:
freeInCell = gpd.sjoin(free, grid)
In [55]:
# Only three cell have free-wifi (specifically, linknyc free wifi)
freeInCell.ID.unique()
Out[55]:
In [79]:
#freeCell = freeInCell.groupby(freeInCell.ID).apply(lambda x: { 'counts': len(x.level), 'mean':mean(x.level), 'median': median(x.level)})
freeCell = freeInCell.groupby(freeInCell.ID).apply(lambda x: median(x.level))
In [84]:
grid2 = grid.copy()
grid2.index = grid2.ID
freeGrid = grid2.join(pd.DataFrame(freeCell, columns=['median level']))
freeGrid.reset_index(drop=True, inplace=True)
freeGrid['median level'].fillna(0, inplace=True)
# freeGrid.sort_values(by='median level')