Bug in gmt_nc.c for grids using 64-bit addressing
|Target version:||Candidate for next bugfix release|
This is a clarification of issue #1019. The key problem is a bug in gmt_nc.c which only appears when the grid exceeds 32-bit addressing. The following example illustrates the problem:
# NC4 format using byte output gmt grdlandmask -r -Rd -Dl -I10s -N0/1 -Gnc4.nc=nb gmt grdimage nc4_chunk.nc=nb -Jx0.0001i -P -Ct.cpt -Anc4.png # Native C using byte C output gmt grdlandmask -r -Rd -Dl -I10s -N0/1 -Gbyte.grd=bb gmt grdimage byte.grd=bb -Jx0.0001i -P -Ct.cpt -Abyte.png
Both commands produce a humongous 129600 by 64800 grid for a total of 7.8 Gb nodes, far beyond the 32-bit 4 Gb limit. To run this, grdimage needs to allocate almost 60 Gb of array space. However, GMT uses 64-bit indexing throughout. Attached are the two png plots produced by grdimage (which uses GDAL to write the PNG images). As you see, the PNG produced by the native byte.grd grid is correct (hence grdlandmask is handling the large grid correctly) while the same command, when written via netCDF, shows typical symptoms of 32/64-bit snafu. I have looked long and hard (and made some changes) to gmt_nc.c but the problem remains. Note: I am assuming the problem lies in GMT since netCDF is used by so many groups for very large data that it is unlikely there is a limitation in netCDF; however, this is an assumption to consider as well.
#1 Updated by Remko almost 2 years ago
I do not see the error either. Do not forget that the bug might be in reading rather than writing.
Have you tried to open nc4_chunk.nc with Panoply?
Also, I do not want to rule out a NetCDF bug here ... it could be in the way NetCDF does chunking.
Have you tried it with chunking switched off?
#3 Updated by Paul almost 2 years ago
- File grab.jpeg added
- Status changed from New to Feedback
This is what the nc4_chunk.nc looks like in Panoply. No chuncking (i.e., classic) makes no difference. So it definitively is in the writing. I also wrote directly via GDAL as netcdf and it looks good. So it has to be our gmt_nc.c code doing something 32-bitish.
- Status changed from Feedback to Resolved
- Assignee set to Paul
- % Done changed from 0 to 100
- Affected version changed from 5.x-svn to 5.4.x-svn
Finally fixed this beast. It was the gmtlib_grd_flip_vertical function that reverses the y-direction of the grid before writing a netCDF file that was the culprit. Fixed in r20031 for both 5.4 branch and trunk.