Skip to content

Add gallery example showing usage of polygon objects from a geopandas.GeoDataFrame (choropleth map) #2796

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 19 commits into from
Nov 10, 2023
Merged
Changes from 11 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
59 changes: 59 additions & 0 deletions examples/gallery/maps/choropleth_map.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,59 @@
"""
Choropleth Map
==============

The :meth:`pygmt.Figure.plot` method allows us to plot geographical data such
as polygons which are stored in a :class:`geopandas.GeoDataFrame` object. Use
:func:`geopandas.read_file` to load data from any supported OGR format such as
a shapefile (.shp), GeoJSON (.geojson), geopackage (.gpkg), etc. You can also
use a full URL pointing to your desired data source. Then, pass the
:class:`geopandas.GeoDataFrame` as an argument to the ``data`` parameter of
:meth:`pygmt.Figure.plot`, and style the geometry using the ``pen`` parameter.
To fill the polygons based on a corresponding column you need to set
``fill="+z"`` as well as select the appropriate column using the ``aspatial``
parameter as shown in the example below.
"""

import geopandas as gpd
import pygmt

# Read polygon data using geopandas
gdf = gpd.read_file("https://geodacenter.github.io/data-and-lab/data/airbnb.zip")

fig = pygmt.Figure()

fig.basemap(
region=gdf.total_bounds[[0, 2, 1, 3]],
projection="M10c",
frame="+tPopulation of Chicago",
)

# The dataset contains different attributes, here we select
# the "population" column to plot.

# First, we define the colormap to fill the polygons based on
# the "population" column.
pygmt.makecpt(
cmap="acton",
series=[gdf["population"].min(), gdf["population"].max(), 10],
continuous=True,
reverse=True,
)

# Next, we plot the polygons and fill them using the defined
# colormap. The target column is defined by the aspatial
# parameter.
fig.plot(
data=gdf[["population", "geometry"]],
pen="0.3p,gray10",
close=True,
Copy link
Member

@yvonnefroehlich yvonnefroehlich Nov 7, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I feel we should explain that the close parameter is used to force closed polygons; maybe already in the introduction text at lines 12-14.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is close required?

Copy link
Member

@yvonnefroehlich yvonnefroehlich Nov 8, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is something I am not 100 % sure about. I orientated myself on the EGU short course example for a choropleth map (https://www.generic-mapping-tools.org/egu22pygmt/ecosystem.html#plotting-geospatial-vector-data-with-geopandas-and-pygmt). And there, the close parameter is set to True to force closed polygons. I already tried it without using this, and the figure looked correct to me.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When plotting polygons without fills, close is required to connect the first and last points of polygons.

When filling polygons, the polygons are always "closed" by default.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Then I would suggest to remove it in the example.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hm. Based on this, it seems to me, we do not need to explicitly set the close parameter to True when creating a choropleth map from polygons provided via a GeoDataFrame.

I just tried to plot the polygons of this GeoDataFrame without any fill. Also this works for me correctly without using close=True:

import geopandas as gpd
import pygmt

gdf = gpd.read_file("https://geodacenter.github.io/data-and-lab/data/airbnb.zip")
fig = pygmt.Figure()

fig.basemap(region=gdf.total_bounds[[0, 2, 1, 3]], projection="M6c", frame="+t ")
fig.plot(data=gdf)

fig.show()

Output figure
geodataframe_withoutfill

Copy link
Member

@seisman seisman Nov 9, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here is an example showing how close works with/without fill.

import pygmt
fig = pygmt.Figure()
data = [[1, 1], [1, 5], [5, 5], [5, 1]]
fig.plot(data=data, projection="X10c/10c", frame=True, region=[0, 10, 0, 10], pen="1p")
fig.shift_origin(xshift="w+1c")
fig.plot(data=data, projection="X10c/10c", frame=True, region=[0, 10, 0, 10], pen="1p", close=True)
fig.shift_origin(xshift="w+1c")
fig.plot(data=data, projection="X10c/10c", frame=True, region=[0, 10, 0, 10], pen="1p", fill="lightblue")
fig.show()

map

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @seisman!

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am wondering whether we should update the code example for the choropleth map of the EGU22 short course and remove the close parameter for consistency, despite the issue that the choropleth map is not displayed correctly. Interestingly, I still do not face this issue when running this Jupyter notebook locally in JupyterLab, neither within a conda environment based on the provided yaml file nor within a conda environment with PyGMT dev version and GMT 6.4.0 installed.

fill="+z",
cmap=True,
aspatial="Z=population",
)


# Add colorbar legend
fig.colorbar(frame="x+lPopulation", position="jML+w6c/0.5c")

fig.show()