How to overlay shapefile data on PyGMT Maps

Utpal Kumar   2 minute read      

Plot the counties of Taiwan using shapefile data on a map created by PyGMT

This time we will see how can we overlay shapefile data on PyGMT map. Here, for example I took the counties data available in shp format from data.gov.tw, and overlay it on the high-resolution map of Taiwan.

Import libraries

We will use the geopandas library to read the shp files.

import pygmt
import os
import geopandas as gpd

Counties shp data

I download the shp data from counties and saved it in countiesData in the working directory. There are multiple files in the countiesData, but we only need the COUNTY_MOI_1090820.shp file. Others are related extention files.

We selected the two counties to highlight on the map - Taipei City, Tainan City.

countiesShp = os.path.join("countiesData","COUNTY_MOI_1090820.shp")

gdf = gpd.read_file(countiesShp)

all_data = []
all_data.append(gdf[gdf["COUNTYENG"]=="Taipei City"])
all_data.append(gdf[gdf["COUNTYENG"]=="Tainan City"])

Plot basemap using PyGMT

Now, we can plot the simple basemap using PyGMT. The benefit of using PyGMT is that we don’t need the coastlines and topographic data separately and the output is high-resolution.

region = [119, 123, 21, 26]

fig = pygmt.Figure()
fig.basemap(region=region, projection="M4i", frame=True)

fig.coast( 
    water='skyblue', 
    shorelines=True)

Overlay the counties

Now, we can overlay the selected counties, fill it with green color, and then put all other counties with the background color (white).

for data_shp in all_data:
    fig.plot(data=data_shp,color="green")
fig.plot(data=countiesShp)

Save the map in raster and vector format

Now, we can save the map in raster and vector format for later use.

fig.savefig('map1.png')
fig.savefig('map1.pdf')
Simple Map of Taiwan with counties
Simple Map of Taiwan with counties

Topographic map

Instead of using simple white as background (which btw looks quite descent to me), we can use the topographic background:

import pygmt
import os
import geopandas as gpd

countiesShp = os.path.join("countiesData","COUNTY_MOI_1090820.shp")

gdf = gpd.read_file(countiesShp)

all_data = []
all_data.append(gdf[gdf["COUNTYENG"]=="Taipei City"])
all_data.append(gdf[gdf["COUNTYENG"]=="Tainan City"])


region = [119, 123, 21, 26]

fig = pygmt.Figure()
fig.basemap(region=region, projection="M4i", frame=True)
fig.grdimage("@srtm_relief_03s", shading=True, cmap='geo')

fig.coast( 
    water='skyblue', 
    shorelines=True)

for data_shp in all_data:
    fig.plot(data=data_shp,color="white", pen=["0.02c", 'white'])
fig.plot(data=countiesShp, pen=["0.02c", 'white'])

fig.savefig('map1.png')
fig.savefig('map1.pdf')
Topographic Map of Taiwan with counties
Topographic Map of Taiwan with counties

Conclusions

We have seen how can we easily add the shapefile data on top of the PyGMT map. We can create simple and topographic map easily using the above code.

Disclaimer of liability

The information provided by the Earth Inversion is made available for educational purposes only.

Whilst we endeavor to keep the information up-to-date and correct. Earth Inversion makes no representations or warranties of any kind, express or implied about the completeness, accuracy, reliability, suitability or availability with respect to the website or the information, products, services or related graphics content on the website for any purpose.

UNDER NO CIRCUMSTANCE SHALL WE HAVE ANY LIABILITY TO YOU FOR ANY LOSS OR DAMAGE OF ANY KIND INCURRED AS A RESULT OF THE USE OF THE SITE OR RELIANCE ON ANY INFORMATION PROVIDED ON THE SITE. ANY RELIANCE YOU PLACED ON SUCH MATERIAL IS THEREFORE STRICTLY AT YOUR OWN RISK.


Leave a comment