I have 2 shapefiles for the UK:
In [3]: # SHAPEFILE 1:
...: # WESTMINISTER PARLIAMENTARY CONSTITUENCY UK SHAPEFILE
...: shapefile1 = "../Westminster_Parliamentary_Constituencies_De
...: cember_2017_UK_BSC_SUPER_SMALL/Westminster_Parliamentary_Constituencies_
...: December_2017_UK_BSC.shp"
In [4]: # SHAPEFILE 2:
...: # LAD19 UK SHAPEFILE
...: shapefile2 = "../03_Maps_March_2020/level3_LAD19_CONTAINS_4_L
...: EVELS_OF_DETAIL/Local_Authority_Districts_December_2019_Boundaries_UK_BU
...: C/Local_Authority_Districts_December_2019_Boundaries_UK_BUC.shp"
In [6]: # LOAD SHAPEFILE 1 INTO GEOPANDAS
...: parl_con = gpd.read_file(shapefile1)
...: parl_con.head()
Out[6]:
FID PCON17CD PCON17NM BNG_E BNG_N LONG LAT Shape__Are Shape__Len geometry
0 11 E14000540 Barking 546099 184533 0.105346 51.5408 5.225347e+07 44697.210277 MULTIPOLYGON (((0.07106 51.53715, 0.07551 51.5...
1 12 E14000541 Barnsley Central 433719 408537 -1.492280 53.5724 1.377661e+08 72932.918783 POLYGON ((-1.42490 53.60448, -1.43298 53.59652...
2 13 E14000542 Barnsley East 439730 404883 -1.401980 53.5391 2.460912e+08 87932.525762 POLYGON ((-1.34873 53.58335, -1.33215 53.56286...
3 14 E14000543 Barrow and Furness 325384 484663 -3.146730 54.2522 8.203002e+08 283121.334647 MULTIPOLYGON (((-3.20064 54.06488, -3.20111 54...
4 15 E14000544 Basildon and Billericay 569070 192467 0.440099 51.6057 1.567962e+08 57385.722178 POLYGON ((0.49457 51.62362, 0.50044 51.61807, ...
In [7]: # SHAPEFILE 1 PROJECTION:
...: parl_con.crs
Out[7]: {'init': 'epsg:4326'}
In [12]: # LOAD SHAPEFILE 2 INTO GEOPANDAS
...: lad19 = gpd.read_file(shapefile2)
...: lad19.head()
Out[12]:
objectid lad19cd lad19nm lad19nmw bng_e bng_n long lat st_areasha st_lengths geometry
0 1 E06000001 Hartlepool None 447160 531474 -1.27018 54.676140 9.684551e+07 50305.325058 POLYGON ((448986.025 536729.674, 453194.600 53...
1 2 E06000002 Middlesbrough None 451141 516887 -1.21099 54.544670 5.290846e+07 34964.406313 POLYGON ((451752.698 520561.900, 452424.399 52...
2 3 E06000003 Redcar and Cleveland None 464361 519597 -1.00608 54.567520 2.486791e+08 83939.752513 POLYGON ((451965.636 521061.756, 454348.400 52...
3 4 E06000004 Stockton-on-Tees None 444940 518183 -1.30664 54.556911 2.071591e+08 87075.860824 POLYGON ((451965.636 521061.756, 451752.698 52...
4 5 E06000005 Darlington None 428029 515648 -1.56835 54.535339 1.988128e+08 91926.839545 POLYGON ((419709.299 515678.298, 419162.998 51...
In [13]: # SHAPEFILE 2 PROJECTION:
...: lad19.crs
Out[13]: {'init': 'epsg:27700'}
With the shapefile using WGS 84 projection, I can successfully plot my choropleth using gv.Polygons:
In [14]: # USE GEOPANDAS DATAFRAME WITH gv.Polygons TO PRODUCE INTERACTIVE CHROPLETH:
...: gv.Polygons(parl_con, vdims='PCON17NM'
...: ).opts(tools=['hover','tap'],
...: width=450, height=600
...: )
Out[14]: :Polygons [Longitude,Latitude] (PCON17NM)\
However if I use the shapefile using OSGB projection then I get an error:
In [15]: # USE GEOPANDAS DATAFRAME WITH gv.Polygons TO PRODUCE INTERACTIVE CHROPLETH:
...: gv.Polygons(lad19, vdims='lad19_name',
...: ).opts(tools=['hover','tap'],
...: width=450, height=600
...: )
DataError: Expected Polygons instance to declare two key dimensions corresponding to the geometry coordinates but 3 dimensions were found which did not refer to any columns.
GeoPandasInterface expects a list of tabular data, for more information on supported datatypes see http://holoviews.org/user_guide/Tabular_Datasets.html
I tried converting the projection used but I just got the same error again when I tried to run gv.Polygons again:
In [16]: lad19.crs
Out[16]: {'init': 'epsg:27700'}
In [17]: lad19.crs = {'init': 'epsg:4326'}
...: lad19.crs
Out[17]: {'init': 'epsg:4326'}
In [19]: # USE GEOPANDAS DATAFRAME WITH gv.Polygons TO PRODUCE INTERACTIVE CHROPLETH:
...: gv.Polygons(lad19, vdims='lad19_name',
...: ).opts(tools=['hover','tap'],
...: width=450, height=600
...: )
DataError: Expected Polygons instance to declare two key dimensions corresponding to the geometry coordinates but 3 dimensions were found which did not refer to any columns.
GeoPandasInterface expects a list of tabular data, for more information on supported datatypes see http://holoviews.org/user_guide/Tabular_Datasets.html
Note that I can successfully plot choropleths for both of these shapefiles using gv.Shape. The only difference using gv.Shape is that with shapefile 1 I don’t need to specify the projection used whereas with shapefile 2 I have to specify crs=ccrs.OSGB()
.
Does anyone know what’s going on here?
Thanks
Shapefile download links:
My issue turned out to be caused by my reprojection step from OSGB to WGS 84.
# THE ORIGINAL PROJECTION ON THE SHAPEFILE
In [16]: lad19.crs
Out[16]: {'init': 'epsg:27700'}
While the result of the following command would suggest that the reprojection step worked
In [17]: lad19.crs = {'init': 'epsg:4326'}
...: lad19.crs
Out[17]: {'init': 'epsg:4326'}
if you look at the geometry attribute you can see that it is still made up of eastings and northings and not longitudes and latitudes as you would expect after reprojecting:
In [8]: lad19["geometry"].head()
Out[8]:
0 POLYGON ((448986.025 536729.674, 453194.600 53...
1 POLYGON ((451752.698 520561.900, 452424.399 52...
2 POLYGON ((451965.636 521061.756, 454348.400 52...
3 POLYGON ((451965.636 521061.756, 451752.698 52...
4 POLYGON ((419709.299 515678.298, 419162.998 51...
Name: geometry, dtype: geometry
The solution was to instead reproject from the original to the desired projection using this method, with the key part being to include inplace=True
:
In [11]: lad19.to_crs({'init': 'epsg:4326'},inplace=True)
...: lad19.crs
Out[11]: {'init': 'epsg:4326'}
The eastings and northings contained in the geometry
column have now been converted to longitudes and latitudes
In [12]: lad19["geometry"].head()
Out[12]:
0 POLYGON ((-1.24098 54.72318, -1.17615 54.69768...
1 POLYGON ((-1.20088 54.57763, -1.19055 54.57496...
2 POLYGON ((-1.19750 54.58210, -1.16017 54.60449...
3 POLYGON ((-1.19750 54.58210, -1.20088 54.57763...
4 POLYGON ((-1.69692 54.53600, -1.70526 54.54916...
Name: geometry, dtype: geometry
and now gv.Polygons can use this shapefile to successfully produce a choropleth map:
In [13]: gv.Polygons(lad19, vdims='lad19nm',
...: ).opts(tools=['hover','tap'],
...: width=450, height=600
...: )
Out[13]: :Polygons [Longitude,Latitude] (lad19nm)