Grdsample problem in GMT 6.0

Q&A
1

Hello. I am a GMT veteran of 25 years, not a newbie.

I have been able to replicate the following problem in GMT 6.0 fresh installs under linux and MacosX.

Say I use grdsample to extract a grid that is higher resolution subset of a larger grid like this:

$gmt grdsample stempmin.6.grd -Gstempmin.grd -R10/12.3/5/16 -I600+n/766+n

where my source grid in this example is a smooth1 degree global geographical grid of surface temperatures with valid floating point values at all nodes:

$ gmt grdinfo stempmin.6.grd
stempmin.6.grd: Title: z
stempmin.6.grd: Command: /home/robd/GMT41/GMT4.1/GMT4.1.3/bin/grdmath /home/robd/AVNGRD/tmp.grd 273.3 SUB = /home/robd/AVNGRD/stempmin.6.grd
stempmin.6.grd: Remark:
stempmin.6.grd: Gridline node registration used [Cartesian grid]
stempmin.6.grd: Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7
stempmin.6.grd: x_min: 0.000 x_max: 360.000 x_inc: 1.000 name: n_columns: 361
stempmin.6.grd: y_min: -90.000 y_max: 90.000 y_inc: 1.000 name: n_rows: 181
stempmin.6.grd: z_min: -55.010 z_max: 36.470 name: z
stempmin.6.grd: scale_factor: 1.000 add_offset: 0.000
stempmin.6.grd: format: classic

The strange problem arises where the maximum longitude value in the small subregion has a fractional part that is less than or equal to 0.5

When that happens, it gives NaN values for all columns >12.0 degrees
$ gmt grd2xyz stempmin.grd | grep NaN | head -1
12.001 16.000 NaN

I can work around the bug by extracting a region where the domain of X has a maximum that has a larger fractional part but this is messy to implement across many scripts:

$ gmt grdsample stempmin.6.grd -Gstempmin.grd -R10/12.6/5/16 -I600+n/766+n
$ gmt grd2xyz stempmin.grd | grep NaN | wc -l
0

If I had been careful about the number of rows and columns I could cut that larger grid in the desired longitude range 10/12.3 to give me correct values at all nodes.

It does not make any difference whether the input/output are defined as global and it doesn’t matter what interpolation parameters I choose: it even gives NaNs using near-neighbor.

It does not matter if the input grid is changed to some other form.

It does not matter if define dx/dy intervals explicitly with -I instead

The bug is not shared with grdtrack, even though grdtrack presumably has much in common with grdsample:

$ echo “12.3 45” | gmt grdtrack -Gstempmin.6.grd
12.300 45.000 8.645
$ echo “12.6 45” | gmt grdtrack -Gstempmin.6.grd
12.600 45.000 10.473

This looks like it should be a trivial bug to fix.

Where I use GMT4.5 there is no problem in grdsample butI have a vague recollection of a similar problem in grdsample around about GMT3.2. Has some ancient buggy code returned?

2

Will fix this later today - can confirm your troubles.

3

Hi again. While I can confirm this problem in 6.0.0 we must have fixed this some time ago since it is not a problem in our curent master repo. We did have a bug in reporting the adjusted w/e/s/n. It said 12 instead of 13; I have fixed that in the master repo as well. Your options are to install from github master or wait patiently for 6.1 later this month.

4

Hi again,

The bug did not go away in 6.1. In mac and linux packages spurious NaNs now appear as latitude rows at the bottom of the same sub-region whereas previously it was the east edge that had these NaNs.

Eg, where I resample a continuous global 0.25 degree grid defined everywhere:

$ gmt grdinfo ht400.6.grd
ht400.6.grd: Title: tmp.grd
ht400.6.grd: Gridline node registration used [Cartesian grid]
ht400.6.grd: Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7
ht400.6.grd: x_min: 0 x_max: 360 x_inc: 0.25 name: x n_columns: 1441
ht400.6.grd: y_min: -90 y_max: 90 y_inc: 0.25 name: y n_rows: 721
ht400.6.grd: z_min: 6148.64990234 z_max: 7670.52978516 name: z
ht400.6.grd: scale_factor: 1 add_offset: 0
ht400.6.grd: format: classic

I get NaNs south of 5.25N:

$ gmt grdsample ht400.6.grd -Ght400.grd -R34.2/43.4308/5.01538/16.7846 -I0.0153846667/0.0153846013 -fg
$  gmt grd2xyz ht400.grd | grep NaN | wc -l
9616
$ gmt grd2xyz ht400.grd | grep NaN | head
34.2	5.24614904575	NaN

It’s the same in my macos version. As in GMT6.0 it remains unrelated to the input grid format or the interpolation scheme. It’s surely caused by a truncation (rounding?) of the input grid region.

5

Sadly, I can confirm your finding. Back to the drawing board.

6

Fix submitted: https://github.com/GenericMappingTools/gmt/pull/3688

7

Now merged into master.

8

Thanks Paul.

Before your fix I already wrote a little tcl script to automatically grdsample over a minimally extended region such that a call to grdcut recovers requested sub-region free from this NaN problem. Is it worth sharing it here? It seems to allows buggy 6.0, 6.1 gmt grdsample lines to be replaced by grdsample.tcl in old scripts with the arguments left unchanged. I ask because I wonder what proportion of GMT users build from src compared to the proportion who wait for package updates? I would guess that a fix like this one would not appear in packages until version 6.2, several months away.

9

Sure, you can post your script here. There is likely to be a 6.1.1 in not too long but timing is as always unclear. Those who build from git just need to pull the latest and type make install so I think that is the simplest way to stay up-to-date and not rely on distros. This is because we do try to fix bugs in the repo as soon as we can and we also continuously add new features to (what will be called) 6.2 eventually, and waiting for a release points void all that good stuff.

10

Save the following as grdsamplefix.tcl and ensure it is an executable on the exe path.
In scripts, replace “gmt grdsample” with “grdsamplefix.tcl” and leave the other args as they are. I have not tested it beyond my own specific problem.

#! /usr/bin/env tclsh

proc Precision {invalue places} {
set formatstring %.$places
append formatstring f
set outvalue [format “$formatstring” $invalue]
return $outvalue
}

proc Nearest {invalue nearest} {
#e.g. round a value to the nearest 50
set tmp [Precision [expr $invalue*(1.0/$nearest)] 0]
return [Precision [expr $tmp*$nearest] 0]
}

set grd [lindex $argv 0]
set otherflags “”
foreach arg [lrange $argv 1 end] {
switch -glob – $arg {
-R* {regsub – {-R} $arg {} region}
-I* {regsub – {-I} $arg {} inc}
-G* {regsub – {-G} $arg {} outgrd}
default {
append otherflags "$arg "
}
}
}
set otherflags [string trim $otherflags]

set Reast [lindex [split $region /] 0]
set Rwest [lindex [split $region /] 1]
set Rsouth [lindex [split $region /] 2]
set Rnorth [lindex [split $region /] 3]

if {[regexp {/} $inc]} {
set xinc [lindex [split $inc /] 0]
set yinc [lindex [split $inc /] end]
} else {
set xinc $inc
set yinc $inc
}

catch {exec gmt grdinfo $grd | grep x_min} grdinfo
set x_min [lindex [split $grdinfo] 2]
set x_max [lindex [split $grdinfo] 4]
set x_inc [lindex [split $grdinfo] 6]
set nx [lindex [split $grdinfo] end]

catch {exec gmt grdinfo $grd | grep y_min} grdinfo
set ny [lindex [split $grdinfo] end]
set y_min [lindex [split $grdinfo] 2]
set y_max [lindex [split $grdinfo] 4]
set y_inc [lindex [split $grdinfo] 6]

set x_min_wide [Nearest $Reast $x_inc]; if {$x_min_wide > $Reast} {set x_min_wide [expr $x_min_wide - $x_inc]}
set xmin $Reast
while {$xmin > $x_min_wide} {
set xmin [expr $xmin - $xinc]
}
set xmin [expr $xmin - $xinc]

set x_max_wide [Nearest $Rwest $x_inc]; if {$x_max_wide < $Rwest} {set x_max_wide [expr $x_max_wide + $x_inc]}
set xmax $Rwest
while {$xmax < $x_max_wide} {
set xmax [expr $xmax + $xinc]
}
set xmax [expr $xmax + $xinc]

set y_min_wide [Nearest $Rsouth $y_inc]; if {$y_min_wide > $Rsouth} {set y_min_wide [expr $y_min_wide - $y_inc]}
set ymin $Rsouth
while {$ymin > $y_min_wide} {
set ymin [expr $ymin - $yinc]
}
set ymin [expr $ymin - $yinc]

set y_max_wide [Nearest $Rnorth $y_inc]; if {$y_max_wide < $Rnorth} {set y_max_wide [expr $y_max_wide + $y_inc]}
set ymax $Rnorth
while {$ymax < $y_max_wide} {
set ymax [expr $ymax + $yinc]
}
set ymax [expr $ymax + $yinc]

set randNum [expr { int(100 * rand()) }]
catch {exec gmt grdsample $grd -R$xmin/$xmax/$ymin/$ymax -I$inc $otherflags -G$randNum.grd} log; puts $log
catch {exec gmt grdcut $randNum.grd -R$region -G$outgrd} log; puts $log
file delete -force $randNum.grd

11

Hi @robd I moved this questions to the Q&A category instead of the Lounge. This is the preferred place to ask questions. :slight_smile: This is mostly to organize things on our side, not something you need to worry about.

12

A potentially related bug in grdsample may not have been fixed. Consider what happens when we grdsample a global grid (0-360) across Greenwich.

First, the input global grid:
localhost:europe robd$ gmt grdinfo tmp1000.6.grd | grep [xy]_min
tmp1000.6.grd: x_min: 0.00000000 x_max: 360.00000000 x_inc: 0.25000000 name: x n_columns: 1441
tmp1000.6.grd: y_min: -90.00000000 y_max: 90.00000000 y_inc: 0.25000000 name: y n_rows: 721

Selection of an area west of Greenwich works:
localhost:europe robd$ gmt grdsample tmp1000.6.grd -R-10/0/0/10 -fg -Gtest.grd
localhost:europe robd$ gmt grdinfo test.grd | grep [xy]_min
test.grd: x_min: -10.00000000 x_max: 0.00000000 x_inc: 0.25000000 (15 min) name: longitude n_columns: 41
test.grd: y_min: 0.00000000 y_max: 10.00000000 y_inc: 0.25000000 (15 min) name: latitude n_rows: 41

Selection of an area east of Greenwich works:
localhost:europe robd$ gmt grdsample tmp1000.6.grd -R0/10/0/10 -fg -Gtest.grd
localhost:europe robd$ gmt grdinfo test.grd | grep [xy]_min
test.grd: x_min: 0.00000000 x_max: 10.00000000 x_inc: 0.25000000 (15 min) name: longitude n_columns: 41
test.grd: y_min: 0.00000000 y_max: 10.00000000 y_inc: 0.25000000 (15 min) name: latitude n_rows: 41

Attempt to select the combined region (fails):
localhost:europe robd$ gmt grdsample tmp1000.6.grd -R-10/10/0/10 -fg -Gtest.grd
localhost:europe robd$ gmt grdinfo test.grd | grep [xy]_min
test.grd: x_min: -10.00000000 x_max: 0.00000000 x_inc: 0.25000000 (15 min) name: longitude n_columns: 41
test.grd: y_min: 0.00000000 y_max: 10.00000000 y_inc: 0.25000000 (15 min) name: latitude n_rows: 41

i.e. it can sample longitudes -10 to 0, and it can sample 0 to +10, but when it samples in the range -10 +10 we only get one side. Consideration of periodic longitude is broken.

13

Code has been updated to handle the additional problem where 6.0, 6.1 grdsample is required to cross 0/360 (if you come across an equivalent problem at -180/180 edit accordingly).

#! /usr/bin/env tclsh

proc Precision {invalue places} {
set formatstring %.$places
append formatstring f
set outvalue [format “$formatstring” $invalue]
return $outvalue
}

proc Nearest {invalue nearest} {
#e.g. round a value to the nearest 50
set tmp [Precision [expr $invalue*(1.0/$nearest)] 0]
return [Precision [expr $tmp*$nearest] 0]
}

catch {exec gmt get FORMAT_FLOAT_OUT} float_format
catch {exec gmt set FORMAT_FLOAT_OUT %.8f}

set grd [lindex $argv 0]
set otherflags “”
foreach arg [lrange $argv 1 end] {
switch -glob – $arg {
-R* {regsub – {-R} $arg {} region}
-I* {regsub – {-I} $arg {} inc}
-G* {regsub – {-G} $arg {} outgrd}
default {
append otherflags "$arg "
}
}
}
set otherflags [string trim $otherflags]

set Rwest [lindex [split $region /] 0]
set Reast [lindex [split $region /] 1]
set Rsouth [lindex [split $region /] 2]
set Rnorth [lindex [split $region /] 3]

puts “$Reast $Rwest”

if {[regexp {/} $inc]} {
set xinc [lindex [split $inc /] 0]
set yinc [lindex [split $inc /] end]
} else {
set xinc $inc
set yinc $inc
}

puts “grdsamplefix.tcl: gmt --version”
catch {exec gmt --version} version
puts “version:$version”
puts “grdsamplefix.tcl: gmt grdinfo $grd -Cn”
catch {exec gmt grdinfo $grd -Cn} grdinfo
puts “grdsamplefix.tcl grdinfo:$grdinfo”

scan $grdinfo “%s %s %s %s %s %s %s %s %s %s %s %s” x_min x_max y_min y_max z0 z1 x_inc y_inc nx ny registration type
puts “grdsamplefix.tcl x_min:$x_min x_max:$x_max x_inc:$x_inc nx:$nx”

set x_min_wide [Nearest $Rwest $x_inc]; if {$x_min_wide > $Rwest} {set x_min_wide [expr $x_min_wide - $x_inc]}
set xmin $Reast
while {$xmin > $x_min_wide} {
set xmin [expr $xmin - $xinc]
}
set xmin [expr $xmin - $xinc]

set x_max_wide [Nearest $Reast $x_inc]; if {$x_max_wide < $Reast} {set x_max_wide [expr $x_max_wide + $x_inc]}
set xmax $Rwest
while {$xmax < $x_max_wide} {
set xmax [expr $xmax + $xinc]
}
set xmax [expr $xmax + $xinc]

set y_min_wide [Nearest $Rsouth $y_inc]; if {$y_min_wide > $Rsouth} {set y_min_wide [expr $y_min_wide - $y_inc]}
set ymin $Rsouth
while {$ymin > $y_min_wide} {
set ymin [expr $ymin - $yinc]
}
set ymin [expr $ymin - $yinc]

set y_max_wide [Nearest $Rnorth $y_inc]; if {$y_max_wide < $Rnorth} {set y_max_wide [expr $y_max_wide + $y_inc]}
set ymax $Rnorth
while {$ymax < $y_max_wide} {
set ymax [expr $ymax + $yinc]
}
set ymax [expr $ymax + $yinc]

set randNum [expr { int(100000 * rand()) }]
if {$Reast > 0 && $Rwest <0} {
#added to handle a second bug where the new region crosses Greenwhich merdian
set x_split1 $Rwest
while {$x_split1 < 0} {
set x_split1 [expr $x_split1 + $xinc]
}
set x_split1 [expr $x_split1 - $xinc]
set x_split2 $Reast
while {$x_split2 > 0} {
set x_split2 [expr $x_split2 - $xinc]
}
set x_split2 [expr $x_split2 + $xinc]
set randNum2 [expr $randNum +1]
set randNum3 [expr $randNum +2]
catch {exec gmt grdsample $grd -R$xmin/$x_split1/$ymin/$ymax -I$inc $otherflags -G$randNum2.grd} log; puts $log
catch {exec gmt grdsample $grd -R$x_split2/$xmax/$ymin/$ymax -I$inc $otherflags -G$randNum3.grd} log; puts $log
catch {exec gmt grdpaste $randNum2.grd $randNum3.grd -G$randNum.grd -fg} log; puts $log
file delete -force $randNum2.grd
file delete -force $randNum3.grd
} else {
puts “grdsamplefix.tcl: gmt grdsample $grd -R$xmin/$xmax/$ymin/$ymax -I$inc $otherflags -G$randNum.grd”
catch {exec gmt grdsample $grd -R$xmin/$xmax/$ymin/$ymax -I$inc $otherflags -G$randNum.grd} log; puts $log
puts “grdsamplefix.tcl: gmt grdcut $randNum.grd -R$region -G$outgrd”
}

catch {exec gmt grdcut $randNum.grd -R$region -G$outgrd} log; puts $log
file delete -force $randNum.grd

catch {exec gmt set FORMAT_FLOAT_OUT $float_format}

14

Hi Paul,

Sorry to report I still (6.1.1) see problems with grdsample crossing zero longitude. The example below would usually be framed as a cut operation (to simplify it I have discarded the accompanying change in grid increments), but it causes no problems in GMT versions before 6.0.

I take a global 1 degree grid in the range 0-360 and show how grdsample onto a domain that includes negative longitudes sends all of these to NaN even though -fg defines the (input) grid as global. -fg is not being entirely ignored because without it the domain of the output would have been longitudes 0 to 50 but neither is it being applied correctly.

localhost:MRFGRD robd$ gmt --version
6.1.1
localhost:MRFGRD robd$ gmt grdinfo ht1000.6.grd
ht1000.6.grd: Title: tmp.grd
ht1000.6.grd: Command: /home/robd/GMT43/bin/grdmath tmp.grd 1 MUL = /home/robd/AVNGRD/ht1000.6.grd
ht1000.6.grd: Remark:
ht1000.6.grd: Gridline node registration used [Cartesian grid]
ht1000.6.grd: Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7
ht1000.6.grd: x_min: 0 x_max: 360 x_inc: 0.25 name: x n_columns: 1441
ht1000.6.grd: y_min: -90 y_max: 90 y_inc: 0.25 name: y n_rows: 721
ht1000.6.grd: z_min: -332.694458008 z_max: 417.657531738 name: z
ht1000.6.grd: scale_factor: 1 add_offset: 0
ht1000.6.grd: format: classic
localhost:MRFGRD robd$ echo "310 50" | gmt grdtrack -Ght1000.6.grd
310	50	-67.4304656982
localhost:MRFGRD robd$ echo "-50 50" | gmt grdtrack -Ght1000.6.grd -fg
-50	50	-67.4304656982
localhost:MRFGRD robd$ gmt grdsample ht1000.6.grd -Gsample.grd -R-50/50/25/50 -fg
localhost:MRFGRD robd$ gmt grd2xyz sample.grd | head -1
-50	50	NaN
15

Can confirm the problem. Can swear it worked when I fixed it but clearly it is still there.