Grdview 3D view

Q&A
1

Dear Kaiō-Samas of the GMT Multiverse,

After reading the grdview man page and several posts dealing with 3D perspective plotting doubts/issues, I am still wondering what I am doing wrong here. My grdview call is:

grdview test.tif -R-84.6333333333/-81.6125/-7.425/-3.39166666667/-6000/0 -JM3i -JZ0.6i -p140/40 -Cturbo -N-6000+ggray -Qi500 -I+a15+ne0.75 > ps

and the output is here:

grdview
grdview841×556 50.3 KB

Is it normal to have that white area between the seafloor relief and the gray sides of the 3D perspective? Please find here the GeoTIFF file in case it is needed. I would appreciate your support to know if this is the normal behavior or it is a known bug in GMT 6.5.0 (yeah, I am not using GMT 6.7.0…).

Many thanks in advance for your nice support as always,

2

Hi @gery,

I can reproduce your issue.

I played around with this. PyGMT (I saw you are using “normal” GMT) allows you to download the built-in GMT remote datasets into an xarray DataArray. This data can then be passed to grdview and plotted in the same way as your GeoTIFF file. Here, I use the Earth relief data.

import pygmt
region_grd = [-84.6333333333, -81.6125, -7.425, -3.39166666667]
region_z = [-84.8, -81.4, -7.6, -3.2, -6000, 0]

frame = ["wSnE", "afg"]
perspective = [130, 30]
shading = "+a15+ne0.75"

# Download elevation grid into an xarray DataArray
grid_ele = pygmt.datasets.load_earth_relief(resolution="30s", region=region_grd)
tif_ele = "test.tif"

for grid in [grid_ele, tif_ele]:

    fig = pygmt.Figure()
    fig.basemap(region=region_z, projection="M10c", frame=frame, perspective=perspective)

    fig.grdview(
        grid=grid,
        cmap="oleron",
        perspective=perspective,
        shading=shading,
        zsize="1.5c",
        surftype="s",
        plane="-7000+ggray",
    )

    fig.colorbar(frame=["x+lelevation", "y+lm"])

    fig.show()

built-in remote dataset

elevation_grid
elevation_grid1920×1175 141 KB

provided geoTiFF file

elevation_geotiff
elevation_geotiff1920×1160 152 KB

Looking at the images above, I feel it’s not supported that the frontal facade follows the exact boundary of your GeoTIFF file, but it is related to a rectangle shape. This probably leads to the strange appearance here.

1 Like
3

Thanks much @yvonnefroehlich!! ich bedanke mich sehr, Sie sind sehr nett!! PyGMT is therefore the way to go, always Python!

All the best,

4

Hi Yvonne. Why is there a white area between the facade and the frame with annotations?

5

The problem seems likely related to -N that gets lost when it starts find NaNs along the grid’s edges. Not sure what to do in these cases. Probably erroring out and refusing to plot the side walls.

6

Because I used a slightly larger region for the figure itself as for the grid:

region_grd = [-84.6333333333, -81.6125, -7.425, -3.39166666667]
region_z = [-84.8, -81.4, -7.6, -3.2, -6000, 0] 

# ...

    fig = pygmt.Figure()
        fig.basemap(region=region_z, projection="M10c", frame=frame, perspective=perspective)

       fig.grdview(
            grid=grid,
# ...
1 Like
7

Thanks - happy to help (:

You can use grdcut to first save the grid for your region and then you should be able to use this file for plotting (e.g., with grdview) in “normal” GMT.

1 Like
8

The problem here is an user error. One cannot ask to plot a grid’s wall when that grid hasn’t it (it has NaNs instead).

9

Totally agree @Joaquim the strange apparence of the plot is due to the missig grid’s wall and that -N or plane does not follow the exact bourdary of the GeoTIFF file.

My idea here was, that one can get an alternative grid but with non-NaN values via the GMT remote datasets and then use this data instead for the map or desired analysis.

1 Like
10

Yes Yvonne, I understood that. But it’s to be clear on the reason of that result. Ofc it would be nice to detect it in advance and refuse to draw the walls but that needs work (and not unnecessarily slow down the rest of the code).