How to color polygons of a geopandas dataframe in pygmt?

Q&A
#1

Hi, I am trying to make a choropleth map from a shapefile using pygmt. I followed one answer on how to plot polygons using pygmt from geodataframe, but when I specify color=gdf.aal_com to color each polygon by column vasues I get an error:

GMTInvalidInput: Can’t use arrays for color if data is matrix or file.

import geopandas as gpd
import pygmt
import os
#%% 
gdf = gpd.read_file('../Albania_loss.shp')
gdf.head()
hex_id pop square_pop NAME_0 aal_com geometry
0 7310 9749.987529 98.742025 Albania 41208.410919 POLYGON ((19.58093 40.88737, 19.52319 40.88737…
1 8112 476.375406 21.826026 Albania 3317.674009 POLYGON ((20.01394 40.13737, 19.95620 40.13737… 
# collect extent region
bound = gdf.total_bounds
region = [bound[0], bound[2], bound[1], bound[3]]

# write to (tmp) file 
# !geopandas (fiona) is not able to overwrite this format
filename = '../albania.gmt'
gdf.to_file(filename, driver='OGR_GMT')

# plot
fig = pygmt.Figure()
fig.basemap(region=region, projection="M20c", frame=True)
fig.plot(data='albania.gmt',color=gdf.aal_com)
fig.show(method='external')

Is there another way to achieve like the following in pygmt?

choro
choro796×389 61.6 KB

#2

Choropleth maps are not easy to do in GMT. I’m able to do them from Julia but I use the -Z option and some helping functions to set the zvals. Not very helpful I know, sorry.

1 Like
#3

Hi @Jamal, you’re in luck, I recently had to do this as well (see Coloring OGR_GMT polygon files based on attribute column). Since you already have an OGR GMT file, you can try this:

fig = pygmt.Figure()
fig.basemap(region=region, projection="M20c", frame=True)
pygmt.makecpt(cmap="batlow", series=[0, 50000, 100], reverse=True, continuous=True)
fig.plot(
    data="albania.gmt",  # input OGR data file
    pen="0.05p,yellow,-",  # set outline colour
    cmap=True,  # use colormap from makecpt
    color="+z",  # color based on the Z column
    close=True,  # force close polygons
    a="Z=aal_com",  # set aspatial column to 'aal_com'
)
fig.colorbar(position="jMR+jMR+w5c+o0.5c/1c+e")
fig.show()

Let us know if that works, might need to tweak a few things depending on your data file. You might find looking at https://docs.generic-mapping-tools.org/latest/plot.html#g helpful.

P.S. We’re looking into integrating pygmt with geopandas at https://github.com/GenericMappingTools/pygmt/issues/608 so if you have any suggestions, feel free to leave a comment there!

1 Like
#4

Thank you @weiji14, this is exactly what I was looking for!
with code [adapted]

output
output465×549 68.5 KB

1 Like