import contextily as ctx
import geopandas as gpd
import pandas as pd
import libpysal as lps
import pysal as ps
import palettable as pltt
import datashader as ds
import datashader.transfer_functions as tf
from datashader.bundling import connect_edges, hammer_bundle
from matplotlib import cm
Data comes from the Office of National Statistics:
from IPython.display import IFrame
url = ('http://geoportal.statistics.gov.uk/'\
'datasets/90e15daaaeef4baa877f4bffe01ebef0_2'\
'?geometry=-3.78%2C52.768%2C-0.27%2C53.346')
print(url)
IFrame(url, 600, 600)
http://geoportal.statistics.gov.uk/datasets/90e15daaaeef4baa877f4bffe01ebef0_2?geometry=-3.78%2C52.768%2C-0.27%2C53.346
shp_link = ('Lower_Layer_Super_Output_Areas_December_2001'\
'_Generalised_Clipped_Boundaries_in_England_and_Wales/'\
'Lower_Layer_Super_Output_Areas_December_2001_'\
'Generalised_Clipped_Boundaries_in_England_and_Wales.shp')
shp_link
'Lower_Layer_Super_Output_Areas_December_2001_Generalised_Clipped_Boundaries_in_England_and_Wales/Lower_Layer_Super_Output_Areas_December_2001_Generalised_Clipped_Boundaries_in_England_and_Wales.shp'
%time db = gpd.read_file(shp_link)\
.set_index('lsoa01cd')\
.to_crs(epsg=3857)
CPU times: user 15.6 s, sys: 138 ms, total: 15.7 s Wall time: 15.8 s
ax = db.plot(alpha=0.5, linewidth=0.3,
color='lime', edgecolor='k',
figsize=(9, 9))
ctx.add_basemap(ax, alpha=0.5,
url=ctx.sources.ST_TONER_BACKGROUND);
db.head()
objectid | lsoa01nm | lsoa01nmw | st_areasha | st_lengths | geometry | |
---|---|---|---|---|---|---|
lsoa01cd | ||||||
E01000001 | 1 | City of London 001A | City of London 001A | 1.298432e+05 | 2370.656583 | POLYGON ((-10549.80993389651 6713902.648684822... |
E01000002 | 2 | City of London 001B | City of London 001B | 2.278982e+05 | 2481.826139 | POLYGON ((-9810.500426142034 6713690.895929562... |
E01000003 | 3 | City of London 001C | City of London 001C | 5.873326e+04 | 1170.012290 | POLYGON ((-10526.31014614473 6714162.541893314... |
E01000004 | 4 | City of London 001D | City of London 001D | 2.292225e+06 | 9764.246314 | POLYGON ((-8735.429215107544 6714066.41961726,... |
E01000005 | 5 | City of London 001E | City of London 001E | 1.891322e+05 | 2202.802920 | POLYGON ((-8612.900660093976 6713330.467610141... |
nodes = pd.DataFrame({'X': db.centroid.x,
'Y': db.centroid.y})\
.reset_index()
nodes.head()
lsoa01cd | X | Y | |
---|---|---|---|
0 | E01000001 | -10773.632739 | 6.713445e+06 |
1 | E01000002 | -10303.097373 | 6.713474e+06 |
2 | E01000003 | -10668.007122 | 6.714103e+06 |
3 | E01000004 | -10361.647705 | 6.712765e+06 |
4 | E01000005 | -8437.087227 | 6.712635e+06 |
canvas = ds.Canvas(plot_height=400, plot_width=400)
agg = canvas.points(nodes, 'X', 'Y')
tf.shade(agg, cmap=cm.viridis)
%time w = lps.weights.Queen.from_dataframe(db)
CPU times: user 43.5 s, sys: 1.82 s, total: 45.4 s Wall time: 47 s
/Users/dani/anaconda/envs/gds/lib/python3.6/site-packages/libpysal/weights/contiguity.py:180: UserWarning: There is one disconnected observation (no neighbors). Island id: 19464 W.__init__(self, neighbors, ids=ids, **kw)
w[0]
{2952: 1.0, 1: 1.0, 2: 1.0, 3: 1.0}
%time edges = w.to_adjlist()
edges.head()
CPU times: user 213 ms, sys: 6.16 ms, total: 219 ms Wall time: 221 ms
focal | neighbor | weight | |
---|---|---|---|
0 | 0 | 2952 | 1.0 |
1 | 0 | 1 | 1.0 |
2 | 0 | 2 | 1.0 |
3 | 0 | 3 | 1.0 |
4 | 1 | 0 | 1.0 |
datashader
%%time
graph0 = connect_edges(nodes, edges,
source='focal', target='neighbor',
x='X', y='Y')
CPU times: user 2.51 s, sys: 72.9 ms, total: 2.59 s Wall time: 2.59 s
canvas = ds.Canvas(x_range=(nodes['X'].min(), nodes['X'].max()),
y_range=(nodes['Y'].min(), nodes['Y'].max()),
plot_height=2000, plot_width=2000)
tf.shade(canvas.line(graph0, 'X','Y', agg=ds.count()))
%%time
graph = hammer_bundle(nodes, edges,
source='focal', target='neighbor',
x='X', y='Y')
CPU times: user 3min 53s, sys: 3min 41s, total: 7min 35s Wall time: 5min 18s
canvas = ds.Canvas(x_range=(nodes['X'].min(), nodes['X'].max()),
y_range=(nodes['Y'].min(), nodes['Y'].max()),
plot_height=3000, plot_width=3000)
tf.shade(canvas.line(graph, 'X','Y', agg=ds.count()))
def w_plot(w, xys, bundled=False,
x_label='X', y_label='Y',
plot_height=400, plot_width=400,
cmap=["lightblue", "darkblue"]):
'''
Generate a bundled plot of the connections stored in a `W` spatial weights matrix
...
Arguments
---------
w : libpysal.W
Spatial weights matrix object from `libpysal`
xys : DataFrame
Table containing the coordinates of the nodes
bundled : Boolean
[Optional. Default=False] Generate bundled or raw representation
x_label : str
[Optional. Default='X'] Name of the column in `xys` containing the X
coordinates of the nodes
y_label : str
[Optional. Default='Y'] Name of the column in `xys` containing the Y
coordinates of the nodes
plot_height : int
[Optional. Default=400] Height of the resulting plot
plot_width : int
[Optional. Default=400] Width of the resulting plot
cmap : list of colors or matplotlib.colors.Colormap, optional
[Optional. Default=["lightblue", "darkblue"]] Matplotlib colormap to use for the graphic
representation
Returns
-------
'''
edges = w.to_adjlist()
graph_builder = connect_edges
if bundled:
graph_builder = hammer_bundle
graph = graph_builder(xys, edges,
source='focal', target='neighbor',
x=x_label, y=y_label)
canvas = ds.Canvas(x_range=(xys[x_label].min(), xys[x_label].max()),
y_range=(xys[y_label].min(), xys[y_label].max()),
plot_height=plot_height, plot_width=plot_width)
return tf.shade(canvas.line(graph, x_label, y_label, agg=ds.count()),
cmap=cmap)
us = gpd.read_file(ps.examples.get_path('NAT.shp'))
us.crs = {'init' :'epsg:4326'}
us = us.to_crs(epsg=3083)
us.plot()
<matplotlib.axes._subplots.AxesSubplot at 0x7fe0d416f1d0>
w = lps.weights.Queen.from_dataframe(us)
xys = pd.DataFrame({'X': us.centroid.x,
'Y': us.centroid.y})
cs = pltt.colorbrewer.sequential.RdPu_3.hex_colors
%time w_plot(w, xys, \
cmap=[cs[0], cs[-1]], \
plot_width=900, plot_height=600)
CPU times: user 656 ms, sys: 31.2 ms, total: 688 ms Wall time: 703 ms
%time w_plot(w, xys, bundled=True, \
cmap=[cs[0], cs[-1]], \
plot_width=3000, plot_height=2000)
CPU times: user 24.9 s, sys: 1.31 s, total: 26.2 s Wall time: 27.6 s
url = 'https://ec.europa.eu/eurostat/cache/GISCO/distribution/v2/nuts/download/ref-nuts-2016-60m.geojson.zip'
'zip://https://ec.europa.eu/eurostat/cache/GISCO/distribution/v2/nuts/download/ref-nuts-2016-60m.geojson.zip'
nuts = gpd.read_file('NUTS_RG_60M_2016_4326_LEVL_3.geojson')
nuts.plot()
<matplotlib.axes._subplots.AxesSubplot at 0x7fe0d435cc88>
Let's keep only those in continental Europe:
eu = nuts.cx[-12:, 30:]\
.reset_index()\
.to_crs(epsg=3857)
eu.plot()
<matplotlib.axes._subplots.AxesSubplot at 0x7fe0cc363438>
w = lps.weights.Queen.from_dataframe(eu)
xys = pd.DataFrame({'X': eu.centroid.x,
'Y': eu.centroid.y})
/home/danie/miniconda3/envs/gds/lib/python3.6/site-packages/libpysal/weights/weights.py:168: UserWarning: There are 28 disconnected observations Island ids: 243, 299, 300, 301, 469, 470, 471, 472, 517, 520, 537, 538, 900, 901, 902, 912, 913, 932, 985, 1247, 1248, 1306, 1345, 1366, 1382, 1383, 1384, 1450 " Island ids: %s" % ', '.join(str(island) for island in self.islands))
cs = pltt.colorbrewer.sequential.Greens_3.hex_colors
%time w_plot(w, xys, \
cmap=[cs[0], cs[-1]], \
plot_width=800, plot_height=900)
CPU times: user 391 ms, sys: 15.6 ms, total: 406 ms Wall time: 441 ms
%%time
out = w_plot(w, xys, bundled=True,\
cmap=[cs[0], cs[-1]], \
plot_width=1900, plot_height=3000)
CPU times: user 10.3 s, sys: 906 ms, total: 11.2 s Wall time: 13.1 s
out
import numpy as np
import matplotlib.pyplot as plt
out = out.data.astype(float)
out[out==0] = np.nan
ax = eu.plot(alpha=0, figsize=(20, 20))
ctx.add_basemap(ax, url=ctx.sources.ST_TONER_BACKGROUND, alpha=0.25)
ext = xys['X'].min(), xys['X'].max(), xys['Y'].min(), xys['Y'].max()
ax.imshow(np.flipud(out), extent=ext, cmap=cm.Greens_r)
ax.set_axis_off()
plt.show()