python-3.xshapefilechoroplethholovizgeoviews

gv.Polygons DataError When Using OSGB Projection


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)\

enter image description here

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:

Shapefile 1:

https://geoportal.statistics.gov.uk/datasets/westminster-parliamentary-constituencies-december-2017-uk-bsc

Shapefile 2:

https://geoportal.statistics.gov.uk/datasets/local-authority-districts-december-2019-boundaries-uk-buc


Solution

  • 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)
    

    enter image description here