Unable to handle GeoDataFrame with mixed Polygons/Multipolygons

[weiji14]

Description of the problem

Originally reported at https://forum.generic-mapping-tools.org/t/not-able-to-draw-closed-polygons-using-pygmt/3663. After some investigation, tt seems that the plot function is unable to plot geopandas.GeoDataFrame objects that contain both Polygon and MultiPolygon shapes. Plotting just Polygon or just MultiPolygon works.

Minimal Complete Verifiable Example

Run on Binder: https://mybinder.org/v2/gh/GenericMappingTools/egu22pygmt/main?urlpath=lab/tree/book/ecosystem.ipynb

# Adapted from https://www.generic-mapping-tools.org/egu22pygmt/ecosystem.html
import pygmt
import geopandas as gpd

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
print(world.head())

This is the first 5 rows of the table

     pop_est      continent                      name iso_a3  gdp_md_est  \
0     920938        Oceania                      Fiji    FJI      8374.0   
1   53950935         Africa                  Tanzania    TZA    150600.0   
2     603253         Africa                 W. Sahara    ESH       906.5   
3   35623680  North America                    Canada    CAN   1674000.0   
4  326625791  North America  United States of America    USA  18560000.0   

                                            geometry  
0  MULTIPOLYGON (((180.00000 -16.06713, 180.00000...  
1  POLYGON ((33.90371 -0.95000, 34.07262 -1.05982...  
2  POLYGON ((-8.66559 27.65643, -8.66512 27.58948...  
3  MULTIPOLYGON (((-122.84000 49.00000, -122.9742...  
4  MULTIPOLYGON (((-122.84000 49.00000, -120.0000...  

When trying to plot the full GeoDataFrame with mixed Polygon/MultiPolygons, it fails

fig = pygmt.Figure()
fig.basemap(region="d", projection="H15c", frame=True)
fig.plot(data=world, pen="1p,black", close=True, color="+z", cmap="bilbao", aspatial="Z=pop_est")
fig.show()

produces

But plotting just MultiPolygon works.

fig = pygmt.Figure()
fig.basemap(region="d", projection="H15c", frame=True)
fig.plot(data=world.iloc[3:5], pen="1p,black", close=True, color="+z", cmap="bilbao", aspatial="Z=pop_est")
fig.show()

produces

I've also documented a workaround using https://geopandas.org/en/v0.12.2/docs/reference/api/geopandas.GeoDataFrame.explode.html to convert all MultiPolygons to Polygons at https://forum.generic-mapping-tools.org/t/not-able-to-draw-closed-polygons-using-pygmt/3663/3, but it's not exactly ideal for some use cases.

Full error message

/srv/conda/envs/notebook/lib/python3.9/site-packages/geopandas/io/file.py:362: FutureWarning: pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.
  pd.Int64Index,
plot [ERROR]: OGR parsing incomplete (is file missing OGR statements?) - abort

System information

PyGMT information:
  version: v0.6.1
System information:
  python: 3.9.15 | packaged by conda-forge | (main, Nov 22 2022, 15:55:03)  [GCC 10.4.0]
  executable: /srv/conda/envs/notebook/bin/python
  machine: Linux-5.10.147+-x86_64-with-glibc2.27
Dependency information:
  numpy: 1.23.5
  pandas: 1.5.2
  xarray: 2022.11.0
  netCDF4: 1.6.0
  packaging: 21.3
  geopandas: 0.10.2
  ghostscript: 9.54.0
  gmt: 6.3.0
GMT library information:
  binary dir: /srv/conda/envs/notebook/bin
  cores: 8
  grid layout: rows
  library path: /srv/conda/envs/notebook/lib/libgmt.so
  padding: 2
  plugin dir: /srv/conda/envs/notebook/lib/gmt/plugins
  share dir: /srv/conda/envs/notebook/share/gmt
  version: 6.3.0

Related

Ideally, there should be a way to read in standard OGR files like GeoJSON directly GenericMappingTools/gmt#4599. Having a mix of Polygons/MultiPolygons is fairly common.

[seisman]

Is it still a bug? The following script works well for me:

import geopandas as gpd
import pygmt

world = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))
world["pop_est"] = world.pop_est * 1e-6

cmap_bounds = pygmt.info(data=world["pop_est"], per_column=True)

fig = pygmt.Figure()
fig.basemap(region="d", projection="H15c", frame=True)
pygmt.makecpt(series=cmap_bounds, cmap="bilbao")
fig.plot(
    data=world,
    pen="1p,black",
    close=True,
    fill="+z",
    cmap=True,
    aspatial="Z=pop_est",
)
fig.colorbar(frame=True, cmap=True)
fig.show()

PyGMT information:
  version: v0.10.1.dev28+g6c933b29a
System information:
  python: 3.11.6 | packaged by conda-forge | (main, Oct  3 2023, 10:40:35) [GCC 12.3.0]
  executable: /opt/miniconda/bin/python3.11
  machine: Linux-6.5.5-200.fc38.x86_64-x86_64-with-glibc2.37
Dependency information:
  numpy: 1.26.0
  pandas: 2.1.1
  xarray: 2023.9.0
  netCDF4: None
  packaging: 23.2
  contextily: None
  geopandas: 0.14.0
  IPython: 8.16.1
  rioxarray: 0.15.0
  ghostscript: 10.02.0
GMT library information:
  binary version: 6.4.0
  cores: 80
  grid layout: rows
  image layout: 
  library path: /opt/miniconda/lib/libgmt.so
  padding: 2
  plugin dir: /opt/miniconda/lib/gmt/plugins
  share dir: /opt/miniconda/share/gmt
  version: 6.4.0

[seisman]

Ping @weiji14

[weiji14]

Hmm, the world map does seem to get plotted properly now:

My version information:

PyGMT information:
  version: v0.10.1.dev75+gd1b01c38.d20231026
System information:
  python: 3.12.0 | packaged by conda-forge | (main, Oct  3 2023, 08:43:22) [GCC 12.3.0]
  executable: /home/user/mambaforge/envs/pygmt/bin/python
  machine: Linux-6.5.0-2-amd64-x86_64-with-glibc2.37
Dependency information:
  numpy: 1.26.0
  pandas: 2.1.2
  xarray: 2023.10.1
  netCDF4: 1.6.5
  packaging: 23.2
  contextily: 1.4.0
  geopandas: 0.14.0
  ipython: None
  rioxarray: 0.15.0
  ghostscript: 10.02.0
GMT library information:
  binary version: 6.4.0
  cores: 16
  grid layout: rows
  image layout: 
  library path: /home/user/mambaforge/envs/pygmt/lib/libgmt.so
  padding: 2
  plugin dir: /home/user/mambaforge/envs/pygmt/lib/gmt/plugins
  share dir: /home/user/mambaforge/envs/pygmt/share/gmt
  version: 6.4.0

Not sure if it's something that has changed with geopandas or between GMT 6.3.0 to 6.4.0 (or PyGMT), but will close this for now.