grdlandmask produces erroneous results for global file at highest resolution
|Status:||In Progress||Start date:||2017-01-05|
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
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).
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.
- % 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.