Not able to draw closed polygons using pygmt

[gilja11]

Hi all,

I am a complete noob regarding pygmt and just started. I would really really use your help. I am trying to understand how to do some basic operations with pygmt.

I am working on my code and am following tutorial on [1]. The main idea is that I have shapefile that I load in the GeoPandas df. Naturally, I have “geometries” column which stores objects such as Polygon and MultyPolygon. I would like to draw these geometries on my basemap, but it is not working. To be specific, when I do

fig.plot(data=world, pen="1p,black", close=True, color="+z", cmap=True, aspatial="Z=pop_est")

I get error "plot [ERROR]: OGR parsing incomplete (is file missing OGR statements?) - abort"

Now, in YT video the author of the tutorial [1] didn’t have any problems, but when you look at [1] you can see also there the same error.
The important snipped of my code is the next one:

fig = pygmt.Figure()
fig.basemap(region=region, projection="M15c", frame=True)
fig.coast(land="black", water="skyblue", borders="1/0.1p,white")

cmap_bounds = pygmt.info(data=shapefile2['Density'], per_column=True)
pygmt.makecpt(cmap="bilbao", series=cmap_bounds)
fig.plot(data = shapefile2, pen="1p,white", close=True, color="+z", aspatial="Z=Density", cmap=True)
fig.show()

If you could help I would be really grateful. Thanks in advance

[1] https://www.generic-mapping-tools.org/egu22pygmt/ecosystem.html

[weiji14]

Hmm, I just noticed that the world map isn’t showing up on the Jupyter Book at Integration with the scientific Python ecosystem 🐍 — Crafting beautiful maps with PyGMT

But it was showing up ok when the recording was made https://youtu.be/72war16Mvxs?t=409 around May 2022. It might be that some updated version of geopandas or some other library has broken this, let us check first.

Edit: Ah yes, I see that the figure has changed at deploy: a4e75825c399a35d8327042c3151d7c0f1a54f91 · GenericMappingTools/egu22pygmt@6b59921 · GitHub. This was how it should have looked like:

Wondering if it’s because geopandas updated their naturalearth_lowres dataset to v5.0.1 at Update naturalearth_lowres dataset to v5.0.1 by rraymondgh · Pull Request #2418 · geopandas/geopandas · GitHub causing some changes to the original dataset we used.

[weiji14]

Ok, it seems like PyGMT/GMT is unable to handle a mix of Polygon/MultiPolygon vector types in a single OGR file, which is the case for the updated naturalearth_lowres v5.0.1 dataset. I’ve submitted a bug report at Unable to handle GeoDataFrame with mixed Polygons/Multipolygons · Issue #2405 · GenericMappingTools/pygmt · GitHub to see if this can be fixed.

As a temporary workaround, you can use geopandas.GeoDataFrame.explode — GeoPandas 0.12.2+0.gefcb367.dirty documentation to convert all the MultiPolygons to Polygons (though this would need to be done carefully). Example code:

import pygmt
import geopandas as gpd

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

# Calculate the populations in millions per capita
world = world[(world.pop_est>0) & (world.name!="Antarctica")]
world['pop_est'] = world.pop_est * 1e-6

# Find the range of data values for creating a colormap
cmap_bounds = pygmt.info(data=world['pop_est'], per_column=True)
cmap_bounds

# Create an instance of the pygmt.Figure class
fig = pygmt.Figure()
# Create a colormap for the figure
pygmt.makecpt(cmap="bilbao", series=cmap_bounds)
# Create a basemap
fig.basemap(region="d", projection="H15c", frame=True)
# Plot the GeoDataFrame
# - Use `close=True` to specify that the polygons should be forced closed
# - Plot the polygon outlines with a 1 point, black pen
# - Set that the color should be based on the `pop_est` using the `color, `cmap`, and `aspatial` parameters
fig.plot(data=world.explode(ignore_index=True), pen="1p,black", close=True, color="+z", cmap=True, aspatial="Z=pop_est")
# Add a colorbar
fig.colorbar(position="JMR", frame='a200+lPopulation (millions)')
# Display the output
fig.show()

should produce: