%matplotlib inline
import pysal
import pandas
import geopandas
import contextily as cx
from h3 import h3
from shapely.geometry import Polygon
H3_ZOOM = 8
url = ("https://github.com/gdsbook/geographic-data-science/"
"raw/master/data/tokyo_clean.csv")
db = pandas.read_csv(url)
ax = db.plot.scatter("longitude", "latitude",
s=2, alpha=0.5,
figsize=(9, 9))
cx.add_basemap(ax, crs={'init': 'epsg:4326'},
url=cx.providers.CartoDB.DarkMatter)
h3
addresses for the extent of the map:¶h3
hexs in the bounding boxe, s, w, n = db.describe()\
.loc[['min', 'max'],
['longitude', 'latitude']]\
.values\
.flatten()
bb_gjson = {'type': 'Polygon',
'coordinates': (((e, s), (w, s),
(w, n), (e, n),
(e, s)),)}
all_hexs = h3.polyfill(bb_gjson,
H3_ZOOM,
geo_json_conformant=True)
len(all_hexs)
1536
polygonise = lambda hex_id: Polygon(
h3.h3_to_geo_boundary(
hex_id, geo_json=True)
)
%time all_polys = geopandas.GeoSeries(list(map(polygonise, all_hexs)), \
index=all_hexs, \
crs={'init': 'epsg:4326'})
CPU times: user 70 ms, sys: 0 ns, total: 70 ms Wall time: 78.8 ms
h3
addresses for areas with pics¶h3fy = lambda r: h3.geo_to_h3(r['latitude'], r['longitude'], H3_ZOOM)
db['h3hex'] = db[['longitude', 'latitude']].apply(h3fy, axis=1)
polys_w_pics = db.groupby('h3hex').size()
polys_w_pics.head()
h3hex 882f5a048dfffff 1 882f5a04b9fffff 1 882f5a04d1fffff 1 882f5a04d5fffff 1 882f5a04d7fffff 1 dtype: int64
polys_db = geopandas.GeoDataFrame({'geometry': all_polys,
'pic_count': 0},
crs=all_polys.crs)
polys_db['pic_count'].update(polys_w_pics)
polys_db.head()
geometry | pic_count | |
---|---|---|
882f5a22d3fffff | POLYGON ((139.634690787546 35.79728757214633, ... | 0 |
882f5a3591fffff | POLYGON ((139.5904945814933 35.764062021256, 1... | 0 |
882f5aaa8dfffff | POLYGON ((139.9094012006458 35.59907375632831,... | 0 |
882f5aae03fffff | POLYGON ((139.7893529828771 35.55657182338908,... | 0 |
882f5aa1bdfffff | POLYGON ((139.8556914797345 35.56658317462952,... | 0 |
ax = polys_db.plot(column='pic_count',
scheme='quantiles', k=3,
alpha=0.25, figsize=(9, 9))
cx.add_basemap(ax, crs=polys_db.crs,
url=cx.providers.CartoDB.DarkMatter);
GeoDatraFrame
%time w_gpd = pysal.lib.weights.Queen.from_dataframe(polys_db)
CPU times: user 250 ms, sys: 0 ns, total: 250 ms Wall time: 258 ms
kRing
def w_from_hids(hids):
shids = set(hids)
neis = {}
for hid in hids:
neis[hid] = list(h3.k_ring(hid, 1).intersection(shids))
w = pysal.lib.weights.W(neis, id_order=hids, ids=hids)
return w
%time w_h3 = w_from_hids(polys_db.index.tolist())
CPU times: user 120 ms, sys: 0 ns, total: 120 ms Wall time: 121 ms