Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Rather than complain, do it? #6588

Closed
joa-quim opened this issue Apr 18, 2022 · 6 comments · Fixed by #6592 or #6615
Closed

Rather than complain, do it? #6588

joa-quim opened this issue Apr 18, 2022 · 6 comments · Fixed by #6592 or #6615
Labels
feature request Request a new feature

Comments

@joa-quim
Copy link
Member

It it detects the cube, it can as well set -Q

grdinfo U.nc
grdinfo [WARNING]: Detected a data cube (U.nc) but -Q not set - skipping
@joa-quim joa-quim added the feature request Request a new feature label Apr 18, 2022
@PaulWessel
Copy link
Member

True, that.

@seisman
Copy link
Member

seisman commented Apr 22, 2022

https://github.com/pydata/xarray-data/raw/master/eraint_uvz.nc

This netcdf file was used in a PyGMT test added by @weiji14. I don't know much about this dataset, but it's likely not a grid cube.

In GMT 6.3.0, GMT reports an error, so it looks reasonable to me:

$ gmt grdinfo -C eraint_uvz.nc -Q
grdinfo [WARNING]: No 3-D array in file eraint_uvz.nc.  Selecting first 3-D slice in the 4-D array z
grdinfo [ERROR]: No 3rd-dimension coordinate vector found in eraint_uvz.nc
grdinfo [ERROR]: gmtapi_import_cube: Unable to examine cube eraint_uvz.nc.
[Session gmt (0)]: Error returned from GMT API: GMT_RUNTIME_ERROR (79)
[Session gmt (0)]: Error returned from GMT API: GMT_RUNTIME_ERROR (79)

In GMT dev (maybe after #6592), it crashes:

$ gmt grdinfo -C eraint_uvz.nc -Vd
....
grdinfo [DEBUG]: Found readable file eraint_uvz.nc
grdinfo [DEBUG]: Object ID 0 : Registered Grid File eraint_uvz.nc as an Input resource with geometry Surface [n_objects = 1]
grdinfo [DEBUG]: gmtapi_begin_io: Input resource access is now enabled [container]
grdinfo [DEBUG]: gmtapi_import_grid: Passed ID = 0 and mode = 1
grdinfo [DEBUG]: Found readable file eraint_uvz.nc
grdinfo [DEBUG]: Calling nc_open on eraint_uvz.nc, ncid = 65536, err = 0
grdinfo [INFORMATION]: No 2-D array in file eraint_uvz.nc.  Selecting first 2-D slice in the 4-D array z
grdinfo [DEBUG]: Calling nc_close on ncid 65536, err = 0
grdinfo [DEBUG]: Calling nc_open on eraint_uvz.nc, ncid = 65536, err = 0
grdinfo [INFORMATION]: No 2-D array in file eraint_uvz.nc.  Selecting first 2-D slice in the 4-D array z
grdinfo [INFORMATION]: No range attribute, guessing registration to be gridline
ERROR: Caught signal number 11 (Segmentation fault: 11) at
0   libgmt.6.4.0.dylib                  0x0000000103d72afc gmtnc_grd_info + 8556
1   ???                                 0x0000600802f15138 0x0 + 105587525374264
Stack backtrace:
0   libgmt.6.4.0.dylib                  0x0000000103df8789 sig_handler_unix + 217
1   libsystem_platform.dylib            0x00007ff81de48e2d _sigtramp + 29
2   libsystem_c.dylib                   0x00007ff85f65e7a0 __global_locale + 0
3   libgmt.6.4.0.dylib                  0x0000000103d43eda gmtlib_read_grd_info + 106
4   libgmt.6.4.0.dylib                  0x0000000103d113e0 gmtapi_import_grid + 5392
5   libgmt.6.4.0.dylib                  0x0000000103d08da8 gmtapi_import_data + 392
6   libgmt.6.4.0.dylib                  0x0000000103cf4dec gmtapi_get_data + 348
7   libgmt.6.4.0.dylib                  0x0000000103cec2ee GMT_Read_Data + 4910
8   libgmt.6.4.0.dylib                  0x0000000103fa225f GMT_grdinfo + 4047
9   libgmt.6.4.0.dylib                  0x0000000103cfd284 GMT_Call_Module + 2180
10  gmt                                 0x00000001036993c8 main + 1816
11  dyld                                0x0000000109f194fe start + 462

@joa-quim
Copy link
Member Author

Apparently not an easy case. When I try to open it with Mirone
image

But GDAL clearly say it's a cube (well, a multi cube file).

gdalinfo eraint_uvz.nc
Driver: netCDF/Network Common Data Format
Files: eraint_uvz.nc
Size is 512, 512
Metadata:
  NC_GLOBAL#Conventions=CF-1.0
  NC_GLOBAL#Info=Monthly ERA-Interim data. Downloaded and edited by [email protected]
Subdatasets:
  SUBDATASET_1_NAME=NETCDF:"eraint_uvz.nc":z
  SUBDATASET_1_DESC=[2x3x241x480] geopotential (16-bit integer)
  SUBDATASET_2_NAME=NETCDF:"eraint_uvz.nc":u
  SUBDATASET_2_DESC=[2x3x241x480] eastward_wind (16-bit integer)
  SUBDATASET_3_NAME=NETCDF:"eraint_uvz.nc":v
  SUBDATASET_3_DESC=[2x3x241x480] northward_wind (16-bit integer)

Though it resists being open with

julia> UVW = gdaltranslate("NETCDF:eraint_uvz.nc:u");
ERROR: BoundsError: attempt to access 480×241×6 Array{Float32, 3} at index [480×241 Matrix{Bool}]
Stacktrace:

and also my HDFexplorer viewer say it's an "not recognizable format"

@PaulWessel
Copy link
Member

I'll try to have a look on the weekend, but must also prepare a seminar for Wednesday in Southampton.

PaulWessel added a commit that referenced this issue Apr 23, 2022
While we correctly found a higher-dimensioned array than 2, we forgot to update ndims once selecting it.
Closes #6588.
PaulWessel added a commit that referenced this issue Apr 23, 2022
* set correct dimension when finding 3- or 4-D array

While we correctly found a higher-dimensioned array than 2, we forgot to update ndims once selecting it.
Closes #6588.

* One more place
@seisman
Copy link
Member

seisman commented May 3, 2022

After PR #6615, the output is:

$ gmt grdinfo eraint_uvz.nc -Cn
-180	179.25	-90	90	66825.5	66825.5	0.75	0.75	480	241	0	1

After PR #6620, the output is:

$ gmt grdinfo eraint_uvz.nc -Cn
grdinfo [WARNING]: No 3-D array in file eraint_uvz.nc.  Selecting first 3-D slice in the 4-D array z
-180	179.25	-90	90	200	850	66825.5	66825.5	0.75	0.75	0	480	241	3	0	1

Can you please verify that the new output is correct?

@seisman seisman reopened this May 3, 2022
@PaulWessel
Copy link
Member

Sure. ncdump -h is all we need:

 ncdump -h eraint_uvz.nc 
netcdf eraint_uvz {
dimensions:
	longitude = 480 ;
	latitude = 241 ;
	level = 3 ;
	month = 2 ;
variables:
	float longitude(longitude) ;
		longitude:_FillValue = NaN ;
		longitude:units = "degrees_east" ;
		longitude:long_name = "longitude" ;
	float latitude(latitude) ;
		latitude:_FillValue = NaN ;
		latitude:units = "degrees_north" ;
		latitude:long_name = "latitude" ;
	int level(level) ;
		level:units = "millibars" ;
		level:long_name = "pressure_level" ;
	short z(month, level, latitude, longitude) ;
		z:number_of_significant_digits = 5 ;
		z:units = "m**2 s**-2" ;
		z:scale_factor = -1.7250274674968 ;
		z:long_name = "Geopotential" ;
		z:add_offset = 66825.5 ;
		z:_FillValue = NaN ;
		z:standard_name = "geopotential" ;
	short u(month, level, latitude, longitude) ;
		u:number_of_significant_digits = 2 ;
		u:units = "m s**-1" ;
		u:scale_factor = -0.00157270493804553 ;
		u:long_name = "U component of wind" ;
		u:add_offset = 26.96875 ;
		u:_FillValue = NaN ;
		u:standard_name = "eastward_wind" ;
	short v(month, level, latitude, longitude) ;
		v:number_of_significant_digits = 2 ;
		v:units = "m s**-1" ;
		v:scale_factor = -0.000477819996337667 ;
		v:long_name = "V component of wind" ;
		v:add_offset = -1.46875 ;
		v:_FillValue = NaN ;
		v:standard_name = "northward_wind" ;
	int month(month) ;

// global attributes:
		:Conventions = "CF-1.0" ;
		:Info = "Monthly ERA-Interim data. Downloaded and edited by [email protected]" ;
}

You can see z, u, and v are all 4-D matrices.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
feature request Request a new feature
Projects
None yet
3 participants