Drapping a grid file over

earthistics · February 5, 2025, 1:46pm

Hello,

I’m trying to drape a grid (.tif) over a topographic grid to visualize it in 2D. However, I’m facing this issue:

→ if I use grdview I get this error message:

grdview [ERROR]: Using this data type (Float64) is not implemented
grdview (gmtapi_import_image): Bad measurement unit.  Choose among c|i|p [/data_wgs84.tif]
[Session pygmt-session (20)]: Error returned from GMT API: GMT_IMAGE_READ_ERROR (22)
[Session pygmt-session (20)]: Error returned from GMT API: GMT_IMAGE_READ_ERROR (22)

→ But if I use grdimage, it successfully generates the plot (using the same .tif file).

Since I can’t obtain any plot with grdview to compare with the obtain with grdimage, I’m wondering what could be the key difference between grdview and grdimage in this case. How could I fix the issue when plotting with grdview?

This is the way I’m trying to use grdview:

region=[lon_min, lon_max, lat_min, lat_max ]
grid_map = pygmt.datasets.load_earth_relief(resolution="01s", region=region)

data_file = "/data_wgs84.tif"

fig = pygmt.Figure()

pygmt.makecpt(cmap='jet', series=[-1.99, 1.99, 0.01])
shade1 = pygmt.grdgradient(grid=grid_map, azimuth="270", normalize="e0.6")
fig.grdimage(grid=grid_map,cmap="gray",shading=shade1)
fig.coast(water='lightblue')
pygmt.makecpt(cmap="gray", series=[0, 3500, 100])
#Draping data over topography
fig.grdview(
    projection="M15c",
    region=region,
    grid=grid_map,
    drapegrid=data_file,
    cmap="jet",
    transparency=25,
    surftype="s",
    frame="a0.2f0.1"
)

pygmt.makecpt(cmap='jet', series=[-1.99, 1.99, 0.01])
fig.colorbar(frame=["xa0.5+lEast", "y+lcm"])

Thanks for your help!

Joaquim · February 5, 2025, 4:14pm

This is more likely a GMT than a PyGMT thing, but not 100% sure. From error message, it’s trying to read the tif image instead of as a grid. So don’t know.

In general GMT doesn’t like Float64’s (for our field they are mostly a waste and too much work to transparently deal with both float types, 32 & 64, in the code), but some attempts are done to be at least be able to ingest them.

Can you post a small example of that .tif file? Use GDAL to cut a chunk of it. You can also try to convert the .tif to 32 bits and see if it works in your case.

earthistics · February 5, 2025, 4:44pm

Thanks for replying. I have converted the .tif file to 32 bits and it still does not work (same error message:

grdview [ERROR]: Using this data type (Float32) is not implemented
grdview (gmtapi_import_image): Bad measurement unit.  Choose among c|i|p

)
Here is a snippet of the file:
chunk_wgs84_f32.tif (78.4 KB)
Thanks for taking a look at this.

Joaquim · February 5, 2025, 7:01pm

Thanks. That shows some GMT bugs and dislike.

dislike. GMT does not like a grid with a .tif extension used as drapegrid. We have to convert it to nc
grdview is bugged when reading a remote @earth_relief file. Must get it first into a file (remote mechanism does a lot of gymns to get the/assemble the data.

That found, this (plain GMR CLI) works.

grdconvert chunk_wgs84_f32.tif -Gchunk_wgs84_f32.grd
grdcut @earth_relief_01s -Rchunk_wgs84_f32.grd -Gt.grd
grd2cpt chunk_wgs84_f32.grd > lixo.cpt
gmt grdview t.grd -JX10 -p210/30 -JZ5 -Qi100 -Ba -Bza -Gchunk_wgs84_f32.grd -Clixo.cpt -png lixo

earthistics · February 5, 2025, 8:20pm

Thanks for having a look at this. Is there any way to convert the .tif file to a .grd one with PyGMT? I didn’t find the PyGMT version of grdconvert.

Joaquim · February 5, 2025, 10:06pm

Having PyGMT means that you have GMT as well but not sure how to get access to it (I’m no PyGMT user). An easy alternative, since you already used it, is to convert the file with GDAL. Note that .grd is just a very old habit (in GMT and other latitudes) to name the extension of netCDF grids. If the file gets named .nc it work equally well.

earthistics · February 6, 2025, 11:45am

Thank you! I converted it to a .nc file and now it works even without cutting the @earth_relief file. However, there’s another issue. I have to handle with NaN values in the drape file and ensure they are transparent, allowing the underlying topography to remain visible. However, it doesn’t work well even when using “c” for the parameter surftype according to the documentation of pygmt.Figure.grdview as you can see:

I’m using grdview like this:

fig.grdview(
    grid=grid_map,
    region=region,
    projection="M12c",  
    frame=["af", "zaf+lElevation (m a.s.l.)"],
    drapegrid=cleaned_grid,
    cmap="disp.cpt",
    zsize="2c",
    surftype="i800",
    shading=True,
    plane="-3000+ggrey",
    interpolation="b",
    perspective=[-200, 35],
    contourpen="0.1p"
)

And this is how it looks without the draped file:

I would appreciate any hint to solve this issue. Thanks!

Joaquim · February 6, 2025, 11:57pm

Not present an immediate answer, and, as I referred PyGMT examples are not for me but search the docs (grdimage is better for that) about transparency modes.