Bug #1019

grdlandmask produces erroneous results for global file at highest resolution

Added by Michael 11 months ago. Updated 8 months ago.

Status:In ProgressStart date:2017-01-05
Priority:NormalDue date:
Assignee:Paul% Done:


Target version:-
Affected version:5.2.x Platform:


Generating a global land/water mask on a 10" grid from the full resolution GSHHG coastlines with grdlandmask gives erroneous results as shown in the attached plots:

Command line used to generate the global mask: grdlandmask -r -Rd -Df -I10s -N0/1 -Gtest.nc?land_flag=nb -V
Note this required about 40 GB of system memory (I used a machine with 60 GB) and ran for several hours.

global.ps shows the problem for a 10 deg x 10 deg box in Indonesia, similar to the one used in the GMT documentation. Note that the mask is ok for part of the region, but completely off for the rest. (I have already verified that this is not a plotting artifact of grdimage by plotting the same region with Panoply.)

For comparison, I did a "regional" run of grdlandmask with same settings otherwise: grdlandmask -r -R125/135/-5/5 -Df -I10s -N0/1 -Gtest1.nc?land_flag=nb -V
This produces correct results, as shown in attached plot regional.ps. I therefore assume that the problem only occurs if the number of grid points to be processed by grdlandmask exceeds a certain threshold.

System information: Redhat Linux 7.3 (kernel 3.10) GMT version 5.2.1
$ uname -a
Linux ip-172-31-39-80.eu-west-1.compute.internal 3.10.0-514.el7.x86_64 #1 SMP Wed Oct 19 11:24:13 EDT 2016 x86_64 x86_64 x86_64 GNU/Linux
$ cat /etc/redhat-release
Red Hat Enterprise Linux Server release 7.3 (Maipo)
$ yum list installed | grep GMT
GMT.x86_64 5.2.1-1.el7 @epel
GMT-common.noarch 5.2.1-1.el7 @epel

global.ps - Results from global grdlandmask run, zooming into Indonesia (582 KB) Michael, 2017-01-05 22:40

regional.ps - Results from regional grdlandmask run for the same region (633 KB) Michael, 2017-01-05 22:41


#1 Updated by Michael 11 months ago

As a workaround, I divided the globe in 4 quadrants and run grdlandmask for each quadrant.

options="-r -Df -I10s -N0/1 -Vl"
gmt grdlandmask -R-180/0/-90/0 -Gsouthwest.nc\?land_flag=nb $options
gmt grdlandmask -R0/180/-90/0 -Gsoutheast.nc\?land_flag=nb $options
gmt grdlandmask -R-180/0/0/90 -Gnorthwest.nc\?land_flag=nb $options
gmt grdlandmask -R0/180/0/90 -Gnortheast.nc\?land_flag=nb $options

I have plotted the resulting land/water masks and found that they were ok, i.e., this avoids the problem reported above which I have seen only when processing global data (-Rd).

#2 Updated by Paul 10 months ago

  • Status changed from New to In Progress
  • Assignee set to Paul

I can confirm the problem and will try to figure out where things go wrong.

#3 Updated by Paul 10 months ago

This is a 32-bit problem, i.e., it fails when the image size exceeds 2^32 pixels, which it does for the global grid. Until we fix it the workaround is to ensure the grid is less than ~4 Gb pixels (as you did). I have fixed some 32-bit bugs in postscriptlight.c and grdimage.c but the problem remains.

#4 Updated by Paul 8 months ago

  • % Done changed from 0 to 30

Some more updates. Some of the GMT grid formats work with 64-bit and the huge grid but presently the default nc4 does not. Still trying to track down that bug. When the formats that do work run grdimage -A to create a png directly then output is correct even at your resolution. However, making the PostScript file may also be compromising since the PostScript language only supports 4-byte signed ints. That does not necessarily mean ghostscript limits us in the same way but it is another possible problem. So tiling seems the smartest option for now.

Also available in: Atom PDF