Bug #755

[GMT5 Ubuntu Linux 14.4 GDAL 2.0.0] Missing GMT driver but can't tell if it's a problem with GMT or GDAL

Added by John over 1 year ago. Updated over 1 year ago.

Status:ClosedStart date:2015-08-21
Priority:NormalDue date:
Assignee:Paul% Done:

100%

Category:-
Target version:Candidate for next bugfix release
Affected version:5.x-svn Platform:

Description

Hi.
In the code below, I'm trying to output a 'classic' GMT grid using the GDAL driver argument that works on my Mac installation. However, I can't remember (or find) an explanation of whether I need to rebuild GMT with GDAL support or GDAL with GMT support (or how to tell which has what support built in or not). Any advice would be welcome while I continue to flail away looking.

J.

+ gmt nearneighbor -N4/1 -V ./Maui_UTM_10m_010.xyz -R738777/764364/2276022/2328753 -I5 -GMaui_UTM_10m_011.nc=gd:GMT -S15
nearneighbor: Grid dimensions are nx = 5118, ny = 10547
nearneighbor: Number of sectors = 4, minimum number of filled sectors = 1
nearneighbor: Processing input table data
nearneighbor: Processed record 13514397
nearneighbor: Gridded row 10547
nearneighbor: GMT_gdalwrite: Output driver GMT not recognized

Qgis-Helly.tiff (118 KB) John, 2015-08-24 13:48

History

#1 Updated by John over 1 year ago

I answered my own question. I rebuilt GDAL with netcdf support and the GMT driver and the NETCDF driver both now show up in the output of:

gdalinfo --formats

J.

#2 Updated by Paul over 1 year ago

  • Status changed from New to Closed
  • Assignee set to Paul
  • % Done changed from 0 to 100

Closing this issue. Not clear to me why you are writing via GDAL when GMT writes netCDF files natively?

#3 Updated by John over 1 year ago

I'm doing it because when I load the output into Qgis, without using the 'classic netcdf' driver, the data appear to be all NaNs. This was the issue that appeared first to me when I was trying to replicate EX24 (I think) and JL explained that since GMT5 the netcdf output is 'chunked' by default. Is there some other way to get around this 'feature'?

#4 Updated by Paul over 1 year ago

Sure. You can turn chunking off and select classic format with the GMT default IO_NC4_CHUNK_SIZE=classic. Or you can write the old GMT netcdf3 format by appending =cf to your filename (c for classic format, f for float). Let me know if the first or the second option works for QGis. I am pretty sure I tried to load a netcdf into Qgis and it worked. What Qgis version do you have?

#5 Updated by John over 1 year ago

Here's the info on my Qgis version. This is from Kyngchaos. The latest version has problems with changing the colormap on OSX so I reverted to this older version which works.

Here's what happened when I tried to resample to 1m with the =c option.

Neptune-Desktop.local:/Volumes/data/mpi-geospatial/src/gmt>./resample.bash
+ INPUT=/Volumes/data/mpi-geospatial/data/level0/Maui_UTM_10m.nc
+ OUTPUT=/Volumes/data/mpi-geospatial/data/level1/Maui-UTM-1m-resampled.nc
+ gmtset PS_MEDIA Custom_7ix5i
+ gmtset PS_PAGE_ORIENTATION portrait
+ INTERVAL=-I1/1
+ grdinfo -I1/1 /Volumes/data/mpi-geospatial/data/level0/Maui_UTM_10m.nc
REGION=-R738772/815549/2276017/2328758
+ PROJECTION=-Ju4/1:3000000
+ TICK=-Ba1.0f1.0g1.0/a1.0f1.0g1.0
+ gmt grdsample /Volumes/data/mpi-geospatial/data/level0/Maui_UTM_10m.nc -G/Volumes/data/mpi-geospatial/data/level1/Maui-UTM-1m-resampled.nc=cf -I1/1
grdsample: GMT_err_trap: -63
grdsample: GMT_err_trap: -63
grdsample (gmt_api.c:2646(GMTAPI_Export_Grid)): NetCDF: Invalid dimension size [/Volumes/data/mpi-geospatial/data/level1/Maui-UTM-1m-resampled.nc=cf]
[Session gmt (0)]: Error returned from GMT API: GMT_GRID_WRITE_ERROR (19)
Neptune-Desktop.local:/Volumes/data/mpi-geospatial/src/gmt>

#6 Updated by Paul over 1 year ago

OK, can you post the output of gmt grdinfo Maui_UTM_10m.nc please. Does smell like a 32-bit issue inside netcdf(?).

#7 Updated by John over 1 year ago

Here you go. Thanks.

grdinfo /Volumes/data/mpi-geospatial/data/level0/Maui_UTM_10m.nc
/Volumes/data/mpi-geospatial/data/level0/Maui_UTM_10m.nc: Title:
/Volumes/data/mpi-geospatial/data/level0/Maui_UTM_10m.nc: Command:
/Volumes/data/mpi-geospatial/data/level0/Maui_UTM_10m.nc: Remark:
/Volumes/data/mpi-geospatial/data/level0/Maui_UTM_10m.nc: Pixel node registration used [Cartesian grid]
/Volumes/data/mpi-geospatial/data/level0/Maui_UTM_10m.nc: Grid file format: ci = GMT netCDF format (32-bit integer, deprecated)
/Volumes/data/mpi-geospatial/data/level0/Maui_UTM_10m.nc: x_min: 738772.797798 x_max: 815548.630544 x_inc: 9.99425055278 name: meters nx: 7682
/Volumes/data/mpi-geospatial/data/level0/Maui_UTM_10m.nc: y_min: 2276017.45467 y_max: 2328757.11484 y_inc: 9.99425055278 name: meters ny: 5277
/Volumes/data/mpi-geospatial/data/level0/Maui_UTM_10m.nc: z_min: 0 z_max: 3052 name: meters
/Volumes/data/mpi-geospatial/data/level0/Maui_UTM_10m.nc: scale_factor: 1 add_offset: 0
/Volumes/data/mpi-geospatial/data/level0/Maui_UTM_10m.nc: format: classic

#8 Updated by Paul over 1 year ago

I think your netcdf library is no good. I recreated your grid size using random numbers:

gmt grdmath -R738772.797798/815548.630544/2276017.45467/2328757.11484 -I9.99425055278 -r 0 3052 RAND = t10.nc
gmt grdinfo t10.nc
t10.nc: Title: Produced by grdmath
t10.nc: Command: grdmath -R738772.797798/815548.630544/2276017.45467/2328757.11484 -I9.99425055278 -r -V 0 3052 RAND = t10.nc
t10.nc: Remark: 
t10.nc: Pixel node registration used [Cartesian grid]
t10.nc: Grid file format: nf = GMT netCDF format (32-bit float), COARDS, CF-1.5
t10.nc: x_min: 738772.797798 x_max: 815548.630544 x_inc: 9.99425055272 name: x nx: 7682
t10.nc: y_min: 2276017.45467 y_max: 2328757.11484 y_inc: 9.99425055334 name: y ny: 5277
t10.nc: z_min: 0.000309821218252 z_max: 3052 name: z
t10.nc: scale_factor: 1 add_offset: 0
t10.nc: format: netCDF-4 chunk_size: 129,129 shuffle: on deflation_level: 3

I then resampled with grdsample:
gmt grdsample t10.nc -Gt.nc -I1/1 -V

which gave
gmt grdinfo t.nc
t.nc: Title: t.nc
t.nc: Command: grdsample t10.nc -Gt.nc -I1/1 -V
t.nc: Remark: 
t.nc: Pixel node registration used [Cartesian grid]
t.nc: Grid file format: nf = GMT netCDF format (32-bit float), COARDS, CF-1.5
t.nc: x_min: 738772.797798 x_max: 815548.630544 x_inc: 0.999997821533 name: x nx: 76776
t.nc: y_min: 2276017.45467 y_max: 2328757.11484 y_inc: 0.999993556504 name: y ny: 52740
t.nc: z_min: -1996.17407227 z_max: 4955.71533203 name: z
t.nc: scale_factor: 1 add_offset: 0
t.nc: format: netCDF-4 chunk_size: 129,129 shuffle: on deflation_level: 3

#9 Updated by John over 1 year ago

Would you pls post your output from nc-config --all? Here's mine.

Neptune-Desktop.local:~/Active/UHMC-ICS418-git>nc-config --all

This netCDF 4.3.2 has been built with the following features:

--cc        → gcc
--cflags → -I/Users/hellyj/src/netcdf-4.3.2/Build/include -I/usr/local/hdf5/include -I/usr/local/zlib-1.2.5/include
--libs → -L/Users/hellyj/src/netcdf-4.3.2/Build/lib -lnetcdf
--has-c++   → no
--cxx →
--has-c++4 → no
--cxx4 →
--fc        → 
--fflags →
--flibs →
--has-f90 → no
--has-dap   > yes
--has-nc2 → yes
--has-nc4 → yes
--has-hdf5 → yes
--has-hdf4 → no
--has-pnetcdf
> no
--prefix    > /Users/hellyj/src/netcdf-4.3.2/Build
--includedir
> /Users/hellyj/src/netcdf-4.3.2/Build/include
--version → netCDF 4.3.2

#10 Updated by Paul over 1 year ago

Seems similar except your version is newer. Perhaps you could run my commands to rule out anything with your maui file?

This netCDF 4.3.3.1 has been built with the following features: 

  --cc        -> /usr/bin/clang
  --cflags    ->  -I/opt/local/include -I/opt/local/include -DNDEBUG
  --libs      -> -L/opt/local/lib -lnetcdf

  --has-c++   -> no
  --cxx       -> 
  --has-c++4  -> no
  --cxx4      -> 

  --fc        -> 
  --fflags    -> 
  --flibs     -> 
  --has-f90   -> no

  --has-dap   -> yes
  --has-nc2   -> yes
  --has-nc4   -> yes
  --has-hdf5  -> yes
  --has-hdf4  -> no
  --has-pnetcdf-> no

  --prefix    -> /opt/local
  --includedir-> /opt/local/include
  --version   -> netCDF 4.3.3.1

#11 Updated by John over 1 year ago

Your version is 4.3.3.1 which is newer. I just built that one and notice that yours used clang and I used gcc. I'll do some more testing and try your code. I probably should rebuild GMT5 with the new netcdf version as well.

#12 Updated by John over 1 year ago

Would you pls repeat your test using the 'classic' output format? I will try also once I rebuild everything.

#13 Updated by Paul over 1 year ago

Done, and I get the same error as you.
grdsample: GMT_err_trap: -63

grdsample: GMT_err_trap: -63
grdsample (gmt_api.c:3270(GMTAPI_Export_Grid)): NetCDF: Invalid dimension size [t.nc=cf]
[Session gmt (0)]: Error returned from GMT API: GMT_GRID_WRITE_ERROR (19)
[Session gmt (0)]: Error returned from GMT API: GMT_GRID_WRITE_ERROR (19)
[Session gmt (0)]: Error returned from GMT API: GMT_GRID_WRITE_ERROR (19)


I will see if I can run in debug and determine what that dimension message means.

#14 Updated by John over 1 year ago

That's weird. I just ran your test and it worked like this.

Neptune-Desktop.local:~/work/scratch>./test.bash
+ gmtset IO_NC4_CHUNK_SIZE=classic
+ gmt grdmath -R738772.797798/815548.630544/2276017.45467/2328757.11484 -I9.99425055278 -r 0 3052 RAND = t10.nc
+ gmt grdsample t10.nc -Gt.nc -I1/1 -V
grdsample: Processing input grid
grdsample: Input grid (738772.797798/815548.630544/2276017.45467/2328757.11484) nx = 7682 ny = 5277 dx = 9.99425055272 dy = 9.99425055334 registration = 1
grdsample: Output grid (738772.797798/815548.630544/2276017.45467/2328757.11484) nx = 76776 ny = 52740 dx = 0.999997821533 dy = 0.999993556504 registration = 1
grdsample: Output grid extrema [-1801/4991.69] exceed extrema of input grid [5.68479e-06/3052]; to clip output use -n...+c

#15 Updated by Paul over 1 year ago

Yes that part work (output is t.nc), but using output as t.nc=bf crashes. Oh, I see you did the IO_NC4 thingy. I will try that too. Meanwhile perhaps you can try -Gt.nc=cf.

#16 Updated by John over 1 year ago

This now works as well. I'll try your test request (-Gt.nc=cf)

./resample.bash
+ INPUT=/Volumes/data/mpi-geospatial/data/level0/Maui_UTM_10m.nc
+ OUTPUT=/Volumes/data/mpi-geospatial/data/level1/Maui-UTM-1m-resampled.nc
+ gmtset IO_NC4_CHUNK_SIZE=classic
+ gmtset PS_MEDIA Custom_7ix5i
+ gmtset PS_PAGE_ORIENTATION portrait
+ INTERVAL=-I1/1
+ grdinfo -I1/1 /Volumes/data/mpi-geospatial/data/level0/Maui_UTM_10m.nc
REGION=-R738772/815549/2276017/2328758
+ PROJECTION=-Ju4/1:3000000
+ TICK=-Ba1.0f1.0g1.0/a1.0f1.0g1.0
+ gmt grdsample /Volumes/data/mpi-geospatial/data/level0/Maui_UTM_10m.nc -G/Volumes/data/mpi-geospatial/data/level1/Maui-UTM-1m-resampled.nc -I1/1 -V
grdsample: Processing input grid
grdsample: Input grid (738772.797798/815548.630544/2276017.45467/2328757.11484) nx = 7682 ny = 5277 dx = 9.99425055278 dy = 9.99425055278 registration = 1
grdsample: Output grid (738772.797798/815548.630544/2276017.45467/2328757.11484) nx = 76776 ny = 52740 dx = 0.999997821538 dy = 0.999993556447 registration = 1
grdsample: Output grid extrema [-7.34316/3052.26] exceed extrema of input grid [0/3052]; to clip output use -n...+c

#17 Updated by John over 1 year ago

That test fails. Seems the gmtset command works but the in-line argument does not.

+ gmt grdmath -R738772.797798/815548.630544/2276017.45467/2328757.11484 -I9.99425055278 -r 0 3052 RAND = t10.nc
+ gmt grdsample t10.nc -Gt.nc=cf -I1/1 -V
grdsample: Processing input grid
grdsample: Input grid (738772.797798/815548.630544/2276017.45467/2328757.11484) nx = 7682 ny = 5277 dx = 9.99425055272 dy = 9.99425055334 registration = 1
grdsample: Output grid (738772.797798/815548.630544/2276017.45467/2328757.11484) nx = 76776 ny = 52740 dx = 0.999997821533 dy = 0.999993556504 registration = 1
grdsample: Output grid extrema [-1898.02/5067.75] exceed extrema of input grid [2.8424e-05/3052]; to clip output use -n...+c
grdsample: GMT_err_trap: -63
grdsample: GMT_err_trap: -63
grdsample (gmt_api.c:2653(GMTAPI_Export_Grid)): NetCDF: Invalid dimension size [t.nc=cf]
[Session gmt (0)]: Error returned from GMT API: GMT_GRID_WRITE_ERROR (19)

#18 Updated by Paul over 1 year ago

Note that the IO_NC4_CHUNK_SIZE=classic writes a netcdf via our new gmt_nc.c code while the =cf writes the old format via gmt_cdf.c. So different code. I will see what I can learn in debug (takes a while to get there).

#19 Updated by Paul over 1 year ago

  • Tracker changed from Support to Bug
  • Status changed from Closed to Resolved
  • Target version set to Candidate for next bugfix release
  • Affected version set to 5.x-svn

So what happens is that we call

        nc_def_dim (ncid, "xysize", header->nm, &xysize_dim);

with header→nm (a size_t) being 4049166240, i.e., nx times ny, and the netCDF function returns error -63: Invalid dimension size. This nm value fits in an unsigned 32-bit integer but NOT in a signed 32-bit integer. Digging further on the netCDF pages, I found this Q/A entry:

Have all netCDF size limits been eliminated?

No, there are still some limits on sizes of netCDF objects, even with the new 64-bit offset format. Each fixed-size variable and the data for one record's worth of a record variable are limited in size to a little less that 4 GiB, which is twice the size limit in versions earlier than netCDF 3.6.

The maximum number of records remains 2^32-1.

The old and deprecated GMT netcdf format stored the 2-D array using a 1-D array and clearly this array cannot exceed 2^31-1. This is the driver you get when you use GDAL's GMT format.
The solution is to add a clear error message in the gmt_cdf.c code so that attempts to use this format for very large grid will result in a clear error message and recommendation to use the new netCDF format. In your case the "classic" format and not by appending =cf or using the GDAL driver. I am adding this check for the =c? mechanism and will look at preventing calls to the gdal driver when the GMT format is chosen and the grid is too large. Will let you know when that is done.

#20 Updated by John over 1 year ago

Awesome. Thanks, Paul.
J.

#21 Updated by Paul over 1 year ago

New error messages have been added to prevent using the old netcdf format (either via =c? or =gd:GMT) when the grid size exceeds 2,147,483,647 nodes (2^31-1), which is the max length of a 1-D netcdf variable. In r14775.

#22 Updated by Paul over 1 year ago

  • Status changed from Resolved to Closed

Closed as fixed.

Also available in: Atom PDF