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!