Problem with adding a shapefile into a gridimage

Q&A PyGMT Q&A
1

I am plotting a shapefile using fig.plot (after reading it with Geopandas) and then adding gridded data using fig.grdimage. Below is my code:

df1=pd.read_csv("/home/jairovs/CMAQ/nox_max_mean.csv")
gdf = gpd.read_file(filename="./ZMVM municipios.shp")
gdf.head()
grid =  pygmt.xyz2grd(x=df1.LON, y=df1.LAT, z=df1.nox_base, spacing=0.01, region=[-99.5,-98.5,19,20])

region=[-99.5,-98.5,19,20]
fig = pygmt.Figure()
fig.coast(land="white", water="cornflowerblue", 
          shorelines="1/0.3p", borders=["1/0.1p", '2/0.9p'],  region=[-99.5,-98.5,19,20], projection="M8i", frame=["afg"])
fig.plot(data=gdf, pen="1p,black")
pygmt.makecpt(cmap="berlin", series=[0, 150, 15])
fig.grdimage(grid=grid, region=[-99.5,-98.5,19,20], projection="M8i", shading= False, cmap=True, transparency= 35)
fig.colorbar(position="JBC+w15h", frame=["af", "y+lNO@-x@- [ppb]"])
fig.show()

If I add only fig.plot and fig.gridimage, the shape file is correctly added to the figure. However, the shapefile disappears when adding fig.coast for showing shorelines, water, etc. I have also tried using fig.basemap instead, but the shapefile does not show. Do you have any idea about what the problem might be and how to solve it?

Thank you very much for creating PyGMT, such a useful tool.

2

Without seeing any of your data or any versions of your map, my first guess would be that fig.coast (with filled land and water) and fig.grdimage should be the first two plotting commands, since these would cover anything plotted previously.

1 Like
3

Hello @JairoVS,

thanks for trying out PyGMT! Great to hear that you find it useful for your work.

Currently, it is difficult to help you, as we don’t know how your data looks like. Please share it with use, to make it easier for people to help you :slightly_smiling_face:. Also showing us the output of gdf.head() as well as the current output figure would be helpful, even it is not looking as you expected.

My initial guess would be, that there is something wrong with the plotting order?

1 Like
4

Hello Yvonne, thank you very much for your suggestions. I am uploading the shapefile I am plotting as fig.plot, the CSV file I plot as grdimage, and a figure of the map I get.

Now, When trying to plot only the shapefile, I get the following error:

psconvert [ERROR]: Unable to decode BoundingBox file (maybe no non-white features were plotted?)

I appreciate any help you can provide on how to solve the issue.

[https://drive.google.com/drive/folders/1ZocFtmbCOEa5D_tirHxLhz6BFetpuNEr?usp=share_link]

5

Your psconvert error means that nothing has been plotted, nothing at all. psconvert cannot detect bounding box of an empty plot.

Your .shp data are projected, https://epsg.io/6362
that is, your .shp data have never been plotted, it cannot be plotted like fig.plot(data=gdf, pen="1p,black"). You need to indicate cartesian projection and a scale for your projected data, as well as the right projection and scale for coast and grdimage. This is doable but not easy or intuitive.

your plot always showed only the result of pygmt.coast and pygmt.grdimage.

you are trying to do approximately what example 28 is doing here https://docs.generic-mapping-tools.org/latest/gallery/ex28.html that is mixing projected data (.shp) together with geographical (coastlines, grid).

the difference between example 28 and your case is that example 28 is simpler as it uses automatic detection of projected frame coordinates using -R@Kilauea.utm.nc (their projected data is grid, not vector .shp like in your case). I think this nice autodetection feature cannot be used with vector data. I tried this before (mixing projected vector data together with geographical coastlines etc) and had to fail back to the manual calculation of the “lower left” and “upper right” corners as in the example 28 for GMT4 https://docs.generic-mapping-tools.org/4/gmt/html/GMT_Docs.html#x1-1630007.28 I don’t know how to do this task in pygmt, that is, to reproject boundaries into geographic coordinates, or to reproject the your .shp data into geographical coordinates.

you can get an overview of your ZMVM.shp like this:

fig.plot(data="./ZMVM.shp", pen="1p,black", projection='X4i',
         region=[2745631.809800,2855436.950900,774927.105200,899488.529100],
         frame="af")

region boundaries are taken from the .shp metadata, by running ogrinfo -al -so ZMVM.shp
note the coordinates on the axes, it is projected (meters), not degrees lon/lat.

you can reproject your whole .shp file to geographical coordinates using ogr2ogr:

ogr2ogr ZMVM_WGS84.shp -t_srs WGS84 ZMVM.shp

and then use ZMVM_WGS84.shp for plotting with pygmt. WGS84 is geographical coordinate system, degrees lon/lat

6

Thank you @mkononets for your detailed answer. It helped me to solve the issue. I tried changing the projection as you suggested to geographical coordinates (WGS84) using geopandas (gpd.to_crs) and then plotting it with pygmt (fig.plot) and it worked! Now I can visualize the three things, fig.coast, fig.grdimage, and fig.plot.I attach a figure of the final map. I appreciate your help!!
NO2_shp|477x500

7

can you please show the working plotting code incl gpd.to_crs conversion using geopandas? like I don’t know much of the python geostuff and examples help learning.

8

Thank you, @mkononets. I am adding the used code to this reply.

#Importing the required libraries
import geopandas as gpd
import pandas as pd
import pygmt

#reading data to be used in the grdimage
df=pd.read_csv("/home/user/CMAQ/nox_max_mean.csv")

#Read the shapefile as a geodataframe using GeoPandas
gdf= gpd.read_file(filename="/home/user/CMAQ/shape_mxc/ZMVM.shp")

#Transform shapefile to a different coordinate system using GeoPandas
gdf2=gdf.to_crs(epsg=4326)

#Creating a grid file from table data
grid =  pygmt.xyz2grd(x=df.LON, y=df.LAT, z=df.nox_base, spacing=0.01, region=[-99.5,-98.5,19,20])

#Making the figure
fig = pygmt.Figure()
fig.coast(land="grey", water="cornflowerblue", 
          shorelines="1/0.3p", borders=["1/0.1p", "2/1p"],  region=[-99.5,-98.5,19,20], projection="M8i", frame=["afg"])
pygmt.makecpt(cmap="berlin", series=[0, 150, 15])
fig.grdimage(grid=grid, region=[-99.5,-98.5,19,20], projection="M8i", shading= False, cmap=True, transparency= 35)
fig.colorbar(position="JBC+w15h", frame=["af", "y+lNO@-x@- [ppb]"])
fig.plot(data=gdf2,region=[-99.5,-98.5,19,20], pen="0.5p,black")
fig.show()