Grdmask and grdsample produce different grids with same region and grid spacing

Q&A
1

Hello, I just upgraded to 6.3.0 (the default supplied with ubuntu 22.04) from 6.0.0 (the default supplied with ubuntu 20.04), and am encountering a new error where grids generated by grdsample and grdmask produce different sized output grids when given the exact same region (-R-119.38/-116.85/33.25/35.08) and grid spacing (-I3c).

First, I create a mask file that I will later use in a grdmath command:

kevin@steel:/tmp/1680195554578-0$ cat mask.xy
-119.38000000000001 34.13
-117.5 33.25
-116.85000000000001 34.19
-118.74999999999999 35.08
-119.38000000000001 34.13
kevin@steel:/tmp/1680195554578-0$ gmt grdmask mask.xy -R-119.38/-116.85/33.25/35.08 -I3c -NNaN/1/1 -Gmask.grd
kevin@steel:/tmp/1680195554578-0$ gmt grdinfo mask.grd 
mask.grd: Title: z
mask.grd: Command: grdmask mask.xy -R-119.38/-116.85/33.25/35.08 -I3c -NNaN/1/1 -Gmask.grd
mask.grd: Remark: 
mask.grd: Gridline node registration used [Cartesian grid]
mask.grd: Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7
mask.grd: x_min: -119.38 x_max: -116.85 x_inc: 0.000833333333333 name: x n_columns: 3037
mask.grd: y_min: 33.25 y_max: 35.08 y_inc: 0.000833333333333 name: y n_rows: 2197
mask.grd: v_min: 1 v_max: 1 name: z
mask.grd: scale_factor: 1 add_offset: 0
mask.grd: format: netCDF-4 chunk_size: 133,130 shuffle: on deflation_level: 3

Then, I use grdsample to resample a lower resolution grid to match:

kevin@steel:/tmp/1680195554578-0$ gmt grdinfo interp_diff.grd 
interp_diff.grd: Title: Produced by grdmath
interp_diff.grd: Command: grdmath base_map.grd interp_resampled.grd ADD = interp_diff.grd
interp_diff.grd: Remark: 
interp_diff.grd: Gridline node registration used [Cartesian grid]
interp_diff.grd: Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7
interp_diff.grd: x_min: -119.38 x_max: -116.85 x_inc: 0.005 name: x n_columns: 507
interp_diff.grd: y_min: 33.25 y_max: 35.08 y_inc: 0.005 name: degree n_rows: 367
interp_diff.grd: v_min: -6.64150381088 v_max: -3.09947729111 name: degree
interp_diff.grd: scale_factor: 1 add_offset: 0
interp_diff.grd: format: netCDF-4 chunk_size: 169,184 shuffle: on deflation_level: 3
kevin@steel:/tmp/1680195554578-0$ gmt grdsample interp_diff.grd -Gtopores_interp_diff.grd -I3c -nl -R-119.38/-116.85/33.25/35.08
kevin@steel:/tmp/1680195554578-0$ gmt grdinfo topores_interp_diff.grd
topores_interp_diff.grd: Title: z
topores_interp_diff.grd: Command: grdsample interp_diff.grd -Gtopores_interp_diff.grd -I3c -nl -R-119.38/-116.85/33.25/35.08
topores_interp_diff.grd: Remark: 
topores_interp_diff.grd: Gridline node registration used [Geographic grid]
topores_interp_diff.grd: Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7
topores_interp_diff.grd: x_min: -119.38 x_max: -116.850833333 x_inc: 0.000833333333333 (3 sec) name: longitude n_columns: 3036
topores_interp_diff.grd: y_min: 33.25 y_max: 35.08 y_inc: 0.000833333333333 (3 sec) name: latitude n_rows: 2197
topores_interp_diff.grd: v_min: -6.64150381088 v_max: -3.09947729111 name: z
topores_interp_diff.grd: scale_factor: 1 add_offset: 0
topores_interp_diff.grd: format: netCDF-4 chunk_size: 132,130 shuffle: on deflation_level: 3

Note that the grid output of the grdsample command has 3036 columns, and the output of the grdmask command has 3307 columns. A potential clue is that the grid created with grdsample says “x_inc: 0.000833333333333 (3 sec)” whereas the grid created with grdmask does not contain the “(3 sec)” annotation.

I have found some other topics with similar issues, e.g. this, but I have node registration. I do notice that the grid generated by grdmask is listed as “Gridline node registration used [Cartesian grid]” and that generated by grdsample is “Gridline node registration used [Geographic grid]”, maybe another clue?

Any ideas what’s going on here and how to avoid this difference?

Thanks!

2

Not sure it is your problem but -I3c ?

3

Looks like -I3c is a the old way to specify arcseconds, according to the current documentation it should be -I3s but it still works the same. Here’s the manpage from back in the GMT 4 days showing c for seconds: https://www.soest.hawaii.edu/gmt/gmt/html/man/grdmask.html

Either way, that doesn’t seem to be the issue. I get the same results if I do -I3s and if I do -I0.000833333333333

4

-I3c and -I3s are the same by backwards compatibility compromise.

5

My guess is a bug. But you could try a few things to see if we can get around it:

  1. Add -fg to grdmask so it flags as a geographic grd
  2. See if entering -R in ddd:mm:ss makes any differences, i.e., -R-119:22:48/…

Since both grids are flagged as gridline registration the difference of one node has to be related to some time round-off going the wrong way.

6

OK, I have a workaround! Adding -fg to those commands didn’t change anything, but it did solve the problem when added to the upstream command to produce the lower resolution grd file used as input to the surface command (interp_diff.grd). Specifying the region in DMS didn’t change anything.

So I have a workaround, but in case this represents a useful test case, I have posted the input data and script to reproduce here: http://opensha.usc.edu/ftp/kmilner/misc/2023_gmt_issue/

Here’s the script that now works:

INC="-I3s"
REG="-R-119.38/-116.85/33.25/35.08"
#REG="-R-119:22:48/-116:51:0/33:15:0/35:4:48"

gmt xyz2grd base_map.xyz -Gbase_map.grd -I0.005 $REG  -D/degree/degree/amp/=/=/=  -:
# do GMT interpolation on the scatter data
gmt surface scatter_diffs.xyz -Ginterpolated.grd -I0.02 $REG  -S20m -T0.0i0.1b -: -h0 2>/dev/null
# resample the interpolated file
gmt grdsample interpolated.grd -Ginterp_resampled.grd -I0.005 $REG -nl

# adding -fg here fixes the issue!
gmt grdmath base_map.grd interp_resampled.grd ADD = interp_diff.grd -fg

# adding -fg to these doesn't change anything
gmt grdmask mask.xy $REG $INC -NNaN/1/1 -Gmask.grd
gmt grdsample interp_diff.grd -Gtopores_interp_diff.grd $INC -nl $REG

gmt grdinfo mask.grd
gmt grdinfo topores_interp_diff.grd

Thanks for pointing me in the right direction!

7

Thanks for posting the data and script. I believe I have fixed the bug - the PR is being reviewed. If interested and able you could build GMT from the PR to see for your self.

8

I am close to what I need in this thread.

I have simple sample gmt bash lat Lon pairs and lines on a globe.

Everything I do gives me a super screen display, I want it to just write to an EPS file

I can’t quite get the correct parameters for gmt image to do it.

Thanks!

I can handle EPS conversion to png or jpg,

Does gmt provide jpg or png as an export option?

9

Looks like you forgot to post what you actually tried and for some reason I am completely unable to guess what you did not write. FYI, GMT modern mode lets you select your output format, e.g., PDF, EPS, JPG etc.

10

Terribly easy to add png and have the png pop out of the bash script!

Added a kill all eog as a crude way to defeat the image popping on screen.

Tnx

11

Please, don’t “pirate” another thread. If you have questions start a new topic … and give enough info ti reproduce the issues.