Figure 6: GeoPandas / spatial data - choropleth map - area with rivers and cities of the states of the USA

[yvonnefroehlich]

This PR adds a JN for GeoPandas - spatial data - GeoDataFrame:

• point geometry - cities
• line geometry - rivers
• polygon geometry - population

Preview:

[gitnotebooks]

Found 1 changed notebook. Review the changes at https://app.gitnotebooks.com/GenericMappingTools/pygmt-paper-figures/pull/6

[seisman]

I'm unsure if we need this example. As I understand it, PyGMT can plot GeoPandas objects, but its functionality for processing them is limited since the aspatial metadata will be lost during processing.

[yvonnefroehlich]

I'm unsure if we need this example. As I understand it, PyGMT can plot GeoPandas objects, but its functionality for processing them is limited since the aspatial metadata will be lost during processing.

I think the idea of having a choropleth map is, to show that there is functionality extending the feature available in GMT.

[seisman]

This example highlights GeoPandas integration in PyGMT. A geopandas.GeoDataFrame may contain any of the standard geometry types: Point, LineString, MultiLineString, Polygon, and MultiPolygon.

The script below plots three different datasets (available from https://www.naturalearthdata.com/) with different geometrie.

• world — country boundaries represented by Polygon and MultiPolygon geometries
• rivers — major rivers represented by LineString geometries
• cities — populated places represented by Point geometries

So, only MultiLineString is not covered.

import pygmt
import geopandas as gpd

world = gpd.read_file("https://naciscdn.org/naturalearth/50m/cultural/ne_50m_admin_0_countries.zip")
rivers = gpd.read_file("https://naciscdn.org/naturalearth/110m/physical/ne_110m_rivers_lake_centerlines.zip")  
cities = gpd.read_file("https://naciscdn.org/naturalearth/110m/cultural/ne_110m_populated_places_simple.zip")

world["POP_EST"] *= 1.0e-5

fig = pygmt.Figure()
fig.basemap(region="=SA", projection="M15c", frame=True)
fig.plot(data=world, pen="0.5p,black", fill="lightgreen")
pygmt.makecpt(cmap="turbo", series=(0, 3000, 100))
fig.plot(
    data=world,
    pen="0.2p,gray10",
    fill="+z",
    cmap=True,
    aspatial="Z=POP_EST",
)
fig.colorbar(frame=True)
fig.plot(data=rivers, pen="1p,blue")
fig.plot(data=cities, style="c0.1c", fill="red", pen="black")
fig.text(x=cities.geometry.x, y=cities.geometry.y, text=cities["name"], font="10p,Helvetica-Bold,black", offset="0.2c/0.2c")
fig.show()

[yvonnefroehlich]

Your example looks quite promising, but I get the following error when running the code:

UnicodeEncodeError: 'charmap' codec can't encode characters in position 0-7: character maps to
encoding with 'cp1252' codec failed

Will look it in detail tomorrow or within the next days.

[seisman]

It's likely because the city names contain non-ASCII characters, which doesn't work well in your console.

This is a likely solution for you:

cities["name"] = cities["name"].apply(lambda x: x.encode("utf-8", errors="ignore").decode("utf-8"))

[yvonnefroehlich]

It's likely because the city names contain non-ASCII characters, which doesn't work well in your console.

This is a likely solution for you:

cities["name"] = cities["name"].apply(lambda x: x.encode("utf-8", errors="ignore").decode("utf-8"))

Oh, thanks! I already tried something similar, but it did solve the problem yet. For now, I focus on the geometry column of the dataframe, then the codes work.

[yvonnefroehlich]

Using region="=SA" shows that country codes are supported. But the overlapping labels for the cities between 5°-25° North do not look nice.

region="SA"	region=[-89, -33, -56.5, 9.7]

[seisman]

What about Africa? Squares are populated cities, while only "worldcity" (city with global significance) are labelled.

import pygmt
import geopandas as gpd

world = gpd.read_file("https://naciscdn.org/naturalearth/50m/cultural/ne_50m_admin_0_countries.zip")
rivers = gpd.read_file("https://naciscdn.org/naturalearth/110m/physical/ne_110m_rivers_lake_centerlines.zip")
cities = gpd.read_file("https://naciscdn.org/naturalearth/110m/cultural/ne_110m_populated_places_simple.zip")
world["POP_EST"] *= 1.0e-6

fig = pygmt.Figure()
fig.basemap(region="=AF", projection="Q15c", frame=True)
pygmt.makecpt(cmap="batlow", series=(0, 270, 100)) 
fig.plot(
    data=world[["POP_EST", "geometry"]],
    pen="1p,gray50",
    fill="+z",
    cmap=True,
    aspatial="Z=POP_EST",
)
fig.colorbar(frame="x+lPopulation (millions)")
fig.plot(data=rivers["geometry"], pen="0.5p,skyblue")
fig.plot(data=cities["geometry"], style="s0.3c", fill="red", pen="1p,black")

mask = cities.worldcity == 1
fig.text(
    x=cities[mask].geometry.x,
    y=cities[mask].geometry.y,
    text=cities[mask]["name"],
    font="10p,Helvetica-Bold,black",
    offset="0c/-0.3c",
    justify="TR",
    fill="white@30",
)
fig.show()
fig.savefig(fname="Fig7_PyGMT_geopandas_ne.pdf")

[yvonnefroehlich]

Really like the right Chicago example, including lake etc.

Thanks! Also like this map. But Natural Earth does not provide such detailed data of Chicago which can be used for the choropleth part. If we want to have the same datasource for all datasets, we probably will go with the USA example showing the area of the different states with rivers and cities. Now need to think about Alaska. It looks a bit strange when just plotting it together with the rest of the USA. Sometimes Alaska is placed seperatly as an inset; maybe we can do this? Edit: See review by Dongdong #6 (comment) .

[yvonnefroehlich]

Depending on how we describe this figure in the manuscript, maybe it makes sense to move the code for insets up to the code part for plotting the polygons of the main map (currently line 18).

[seisman]

That would change the region and projection settings and we have to set region and projection again when plotting rivers and cities.

[yvonnefroehlich]

Yes, that's true. We do not use Figure.inset here.

[seisman]

Maybe use a different CPT or a different range, to avoid polygons that are too light or too dark?

[yvonnefroehlich]

Actually not so easy, as Alaska is so much larger and Hawaii so much smaller than the other states.

hawaii	nuuk	imola	bamako

acton	batlow	navia	bilbao (currently)

[seisman]

Maybe "hawaii" is better?

[yvonnefroehlich]

Yes, from the CPTs I tested I think "hawaii" is the best (even I am not really a fan of the colors in this colormap itself).
I have chaged the colormap in commit e1121f2. Also played around with a logarithmic scaling, but did not improve the image.

[yvonnefroehlich]

Just want to let you know, I am considering adding the code for the Chicago map to my personal map repo (please see yvonnefroehlich/gmt-pygmt-plotting#131). As this example relies on a mixture of data sources and uses the geodatasets, I think we will not consider it for the paper or for a tutorial on GeoPandas geometries. So hope this is fine for everybody.

[seisman]

This PR looks good to me.

[seisman]

@yvonnefroehlich If you feel this example is good, please copy the codes to Notebook, remove the Python script, output image, then we should be able to merge this PR.