One of the grdcut bugs

[joa-quim]

The bug shown in

grdcut @earth_relief_03s_g -R-10/-9.8/4.9/5.0 -Glixo.grd

Was introduced in the commit

4/7/2023 SHA-1: bf4918172f8f627963180898334740887e67f42a

• allow slop in wesn (#7356)

But note that that commit only made appear the peak at the SW corner, the one at the NW is cause by a much older code (don't know more).

[joa-quim]

The other bug is (it crashes).

gmt grdcut @earth_relief_03s_g -R-11/-9/4/6.0 -Glixo.grd
grdblend [ERROR]: Failed to remove c:\TMP/grdblend_resampled_a28116.nc! [remove error: Permission denied]
grdblend [ERROR]: Failed to delete file c:\TMP/grdblend_resampled_a28116.nc

Hmm, doesn't crash in the debugger but the result is this.

[PaulWessel]

Strange, I got no error and this plot:

gmt grdcut @earth_relief_03s_g -R-11/-9/4/6.0 -Glixo.grd
gmt grdimage lixo.grd -B -I+d -png a

But for some reason on Windows the removal of the temp file is not allowed/

[joa-quim]

Have to leave right now, but the permission to remove is not the problem . That happens because the struct has twice the same file name and the error message come from the the first time it tries to remove it but the it's still open in grdcut.

[seisman]

Maybe it's also related to https://github.com/GenericMappingTools/pygmt/issues/2511?

[joa-quim]

I don't get the second error on WSL either.

[joa-quim]

And there is also #7012

[PaulWessel]

First bug issue. So grdcut calls grdblend to make a blend of 15s pixel and 3s gridline. The z-range of the entire 1x1 degree 3s tile is -1 to 33 m and inside the selected area it is all zeros. The 15s grid in this area has z = -2048.5 to -940. So, what happens is that grdblend runs grdsample on the 15s subset to bring it onto the 3s nodes, and then we blend, which of course is adding 0 to the 15s value and taking the average so you think it would be close to -1000. Yet the lixo grid values are about or matches matches the 3s grid.

I attach a script that shows a bunch of plots and has those two peaks in the corners.

#!/bin/bash
# Cut&blend the 3s and 15s subsets
gmt grdcut @earth_relief_03s_g -R-10/-9.8/4.9/5.0 -Glixo03s.grd
gmt grdcut @earth_relief_15s_p -R-10/-9.8/4.9/5.0 -Glixo15s.grd
# Cut the raw grids
gmt grdcut ~/.gmt/server/earth/earth_relief/earth_relief_03s_g/N04W010.earth_relief_03s_g.nc -R-10/-9.8/4.9/5.0 -Glixo03s_raw.grd
gmt grdcut ~/.gmt/server/earth/earth_relief/earth_relief_15s_p/N00W010.earth_relief_15s_p.nc -R-10/-9.8/4.9/5.0 -Glixo15s_raw.grd
# The 3s grid blend of 3 and 15s data
gmt begin blend
    gmt grdview lixo03s.grd -R-10/-9.75/4.85/5/-2200/-900 -Bxaf -Byaf -Bzaf -BWSZ+t"Blend" -JZ2i -p135/35 --MAP_TITLE_OFFSET=32p
gmt end show
# Just the 3s grid
gmt begin map3s
    gmt grdview lixo03s_raw.grd -R-10/-9.75/4.85/5/-2200/-900 -Baf -Bzaf -BWSZ+t"3s only" -JZ2i -p135/35 --MAP_TITLE_OFFSET=112p
gmt end show
# Just the 15s grid
gmt begin map15s
    gmt grdview lixo15s_raw.grd -R-10/-9.75/4.85/5/-2200/-900 -Baf -Bzaf -BWSZ+t"15s only" -JZ2i -p135/35 --MAP_TITLE_OFFSET=0p
gmt end show
# Plot the two node sets for the full requested area
gmt begin all
    gmt grdmath -R-10/-9.8/4.9/5.0 -I15s -rp X = t.nc
    gmt grd2xyz t.nc | gmt plot -R-10/-9.75/4.85/5 -Sc2p -Glightblue -B
    gmt grd2xyz lixo03s.grd | gmt plot -Sc1p -Gred
    gmt grdinfo -Ib t.nc | gmt plot -Wfaint
gmt end show
# Plot the two node sets for a zoom in
gmt begin zoom
    gmt grd2xyz t.nc | gmt plot -R-10:0:05/-9:59/4:53:50/4:54:30 -Sc5p -Gblue -Ba9sf3sg3s
    gmt grdinfo -Ib t.nc | gmt plot -W1p
    gmt grd2xyz lixo03s.grd | gmt plot -Sc3p -Gred
gmt end show

[PaulWessel]

So not hard to understand the spike at the SW corner now. See the original node layout with blue being the 15s ocean and red the 3s SRTM. Since the 15s needs to be resampled at the red nodes there are lots of nodes "outside" the 15s data domain and here we get whatever the BCR interpolation feels like. Since there is a huge difference in the values of 15s and 3s in this entire region we get interesting results. I think the only way this would work well is if the 15s region extended outside the 3s region so that we had blue nodes outside to guide the interpolation. So maybe this is a bug or maybe it is a site combination of two largely incompatible data sets sampled hear their edges.

[PaulWessel]

Note if we extend the region to start at -10.25 we still get those two spikes, now in the middle of the grid/map.

[joa-quim]

See my first post. The SW spike started to appear after commit 4/7/2023 SHA-1: https://github.com/GenericMappingTools/gmt/commit/bf4918172f8f627963180898334740887e67f42a which was made to fix the differences in grid sizes with two different modules.

[PaulWessel]

A bit more forensic work. In our case (given the -R), we have data beyond the S, E, and N bounds. Hence, the pad outside those 3 sides are filled with data values from the tile. On the W side we are at the tile boundary (-10) and hence there are no further data values n that tile we can use. Thus, after we read in the 15s sub grid it looks like this, with the fat lines indicating the grid boundary and things outside are the pad.

D means we have a data value, and BC means we need to fill in the pad using the BCs. In our particular case:

• For BC1 nodes in the first boundary column we use the zero normal curvature to solve for those nodes based on the two on the inside (all in x-direction)
• For BC2 nodes we set dLaplacian/dn = 0 and solve using the data and the now set BC1 nodes.
• That step sets the two blue nodes as well, of course.
• To set the green corner nodes we use the d2z/dxdy = 0 setting. Note that this will require the values at the blue nodes.
• There are no BCs controlling what the red and orange cells should be in absence of actual data. Hence these should be set to NaN so that the bcr convolution will skip them.

I think these are the bugs:

1. We compute the green corners before the blue nodes have been assigned, and hence those are initially zero.
2. There are no assignments to NaN, meaning the red and orange nodes default to zero.

Convolving over the top left 4x4 patch thus involves a lot of zeros, pulling the resulting extrapolation up into the spikes.

[seisman]

Ping @PaulWessel to revisit the issue and https://github.com/GenericMappingTools/pygmt/issues/2511. Hopefully they can be fixed before releasing GMT 6.5.

[PaulWessel]

Will try