#Brexit
¶Sources:
%matplotlib inline
import matplotlib.pyplot as plt
import geopandas as gpd
import pandas as pd
import pysal as ps
import urllib2
The official referendum results are available from the Electoral Commission:
ref_link = "http://www.electoralcommission.org.uk/__data/assets/file/0014/212135/EU-referendum-result-data.csv"
This is a csv
that may be read directly from the EC server:
ref = pd.read_csv(ref_link, index_col='Area_Code')
Results are available at the Local Authority level, which we can link to the ONS official boundaries. The GeoJSON
needs to be downloaded locally to read it:
lad_link = 'http://opengeography.ons.opendata.arcgis.com/datasets/943be9957ccf419399bfafc2ad92bc0d_0.geojson'
# Download the file
buf = urllib2.urlopen(lad_link)
lads = buf.read()
tmp = open('lads.geojson', 'w')
tmp.write(lads)
tmp.close()
Now we can read it up:
lads = gpd.read_file('lads.geojson', driver='GeoJSON').set_index('LAD14CD')
And plot it:
lads.plot(linewidth=0);
And now we can simply link the two sides, referendum and spatial data:
refl = gpd.GeoDataFrame(lads.join(ref), crs=lads.crs)\
.to_crs(epsg=4277)
# Write out to a temporaty shapefile
refl.to_file('tmp')
# Create weights
wq = ps.queen_from_shapefile('tmp/tmp.shp')
# This might not work on Windows
! rm -r tmp
WARNING: there are 6 disconnected observations Island ids: [45, 51, 331, 339, 342, 358]
Because of islands, we create a Knn-1 matrix as well:
pts = refl.centroid
xys = pd.np.array([(pt.x, pt.y) for pt in pts])
knn1 = ps.knnW_from_array(xys, k=1)
And now we join both using a union
:
w = ps.weights.Wsets.w_union(wq, knn1)
w.transform = 'R'
lisa = ps.Moran_Local(refl['Pct_Leave'].values, w, permutations=999)
Now let us pull out areas with statistically significant spatial clustering (at the 5% level):
refl['significant'] = lisa.p_sim < 0.05
refl['quadrant'] = lisa.q
At this point, we are ready to visualize the output!
# Setup the figure and axis
f, ax = plt.subplots(1, figsize=(18, 18))
# Plot HH clusters
hh = refl.loc[(refl['quadrant']==1) & (refl['significant']==True), 'geometry']
hh.plot(ax=ax, color='red', linewidth=0.1, edgecolor='w')
# Plot LL clusters
hh = refl.loc[(refl['quadrant']==3) & (refl['significant']==True), 'geometry']
hh.plot(ax=ax, color='blue', linewidth=0.1, edgecolor='w')
# Plot LH clusters
hh = refl.loc[(refl['quadrant']==2) & (refl['significant']==True), 'geometry']
hh.plot(ax=ax, color='#83cef4', linewidth=0.1, edgecolor='w')
# Plot HL clusters
hh = refl.loc[(refl['quadrant']==4) & (refl['significant']==True), 'geometry']
hh.plot(ax=ax, color='#e59696', linewidth=0.1, edgecolor='w')
# Non-significant
ns = refl.loc[refl['significant']!=True, 'geometry']
ns.plot(ax=ax, color='0.75', linewidth=0.1, edgecolor='w')
# Style and draw
f.suptitle('LISA of % for Leave', size=30)
f.set_facecolor('1')
ax.set_axis_off()
plt.axis('equal')
plt.show()