search

GenericMappingTools / gmt

The Generic Mapping Tools
https://www.generic-mapping-tools.org
Other
847 stars 355 forks source link

gmt grdinfo not reporting z_min and z_max correctly for some netCDF files. #8586

Closed hemantbist closed 2 weeks ago

hemantbist commented 2 weeks ago

grdinfo reports z_min and zmax incorrectly as 0. gdalinfo -stats correctly shows the minimum and maximum elevation values (STATISTIC* ). This happens with only some netCDF files ( For most other netCDF / TIF grid files , it work fine. ) I have tried it with 6.0.0 and 6.4.0 version of gmt tools.

Sample netCDF file is attached as a.nc.zip.

Relevant lines from output (The commands and their output also pasted at end) ---- grdinfo unable to report the elevation stats a.nc: z_min: 0 z_max: 0 name: GDAL Band Number 1 ---- gdalinfo correctly reports stats STATISTICS_MAXIMUM=1058.4227294922 STATISTICS_MEAN=1037.8334030559 STATISTICS_MINIMUM=1010.0650634766 STATISTICS_STDDEV=11.264772263914

$gmt grdinfo a.nc a.nc: Title: GDAL Band Number 1 a.nc: Command: Fri Sep 27 09:18:14 2024: GDAL CreateCopy( /tmp/a.nc, ... ) a.nc: Remark: a.nc: Gridline node registration used [Cartesian grid] a.nc: Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7 a.nc: x_min: 710976.648251 x_max: 714542.280646 x_inc: 0.999616595327 name: n_columns: 3568 a.nc: y_min: 4287917.95093 y_max: 4290524.95101 y_inc: 0.999616593282 name: n_rows: 2609 a.nc: z_min: 0 z_max: 0 name: GDAL Band Number 1 a.nc: scale_factor: 1 add_offset: 0 a.nc: format: classic $gdalinfo -stats a.nc Driver: netCDF/Network Common Data Format Files: a.nc a.nc.aux.xml Size is 3568, 2609 Coordinate System is: PROJCS["NAD83(2011) / UTM zone 10N", GEOGCS["NAD83(2011)", DATUM["NAD83_National_Spatial_Reference_System_2011", SPHEROID["GRS 1980",6378137,298.2572221010042, AUTHORITY["EPSG","7019"]], AUTHORITY["EPSG","1116"]], PRIMEM["Greenwich",0], UNIT["degree",0.0174532925199433], AUTHORITY["EPSG","6318"]], PROJECTION["Transverse_Mercator"], PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",-123], PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",500000], PARAMETER["false_northing",0], UNIT["metre",1, AUTHORITY["EPSG","9001"]], AUTHORITY["EPSG","6339"]] Origin = (710976.148442213074304,4290525.450813457369804) Pixel Size = (0.999616595327095,-0.999616593282464) Metadata: Band1#BandName=julall_nad_wg Band1#grid_mapping=transverse_mercator Band1#long_name=GDAL Band Number 1 Band1#RepresentationType=ATHEMATIC Band1#_FillValue=-3.4028235e+38 NC_GLOBAL#Conventions=CF-1.5 NC_GLOBAL#GDAL=GDAL 2.4.1, released 2019/03/15 NC_GLOBAL#GDAL_AREA_OR_POINT=Area NC_GLOBAL#GDAL_DataType=Generic NC_GLOBAL#history=Fri Sep 27 09:18:14 2024: GDAL CreateCopy( /tmp/a.nc, ... ) transverse_mercator#false_easting=500000 transverse_mercator#false_northing=0 transverse_mercator#GeoTransform=710976.1484422131 0.9996165953270648 0 4290525.450813457 0 -0.9996165932824028 transverse_mercator#grid_mapping_name=transverse_mercator transverse_mercator#inverse_flattening=298.2572221010042 transverse_mercator#latitude_of_projection_origin=0 transverse_mercator#longitude_of_central_meridian=-123 transverse_mercator#longitude_of_prime_meridian=0 transverse_mercator#long_name=CRS definition transverse_mercator#scale_factor_at_central_meridian=0.9996 transverse_mercator#semi_major_axis=6378137 transverse_mercator#spatial_ref=PROJCS["NAD83(2011) / UTM zone 10N",GEOGCS["NAD83(2011)",DATUM["NAD83_National_Spatial_Reference_System_2011",SPHEROID["GRS 1980",6378137,298.2572221010042,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","1116"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","6318"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","6339"]] x#long_name=x coordinate of projection x#standard_name=projection_x_coordinate x#units=m y#long_name=y coordinate of projection y#standard_name=projection_y_coordinate y#units=m Corner Coordinates: Upper Left ( 710976.148, 4290525.451) (120d34'21.48"W, 38d44'17.62"N) Lower Left ( 710976.148, 4287917.451) (120d34'24.34"W, 38d42'53.08"N) Upper Right ( 714542.780, 4290525.451) (120d31'53.87"W, 38d44'14.52"N) Lower Right ( 714542.780, 4287917.451) (120d31'56.78"W, 38d42'49.99"N) Center ( 712759.464, 4289221.451) (120d33' 9.12"W, 38d43'33.81"N) Band 1 Block=3568x1 Type=Float32, ColorInterp=Undefined Min=1010.065 Max=1058.423 Minimum=1010.065, Maximum=1058.423, Mean=1037.833, StdDev=11.265 NoData Value=-3.4028234663852886e+38 Metadata: BandName=julall_nad_wg grid_mapping=transverse_mercator long_name=GDAL Band Number 1 NETCDF_VARNAME=Band1 RepresentationType=ATHEMATIC STATISTICS_MAXIMUM=1058.4227294922 STATISTICS_MEAN=1037.8334030559 STATISTICS_MINIMUM=1010.0650634766 STATISTICS_STDDEV=11.264772263914 _FillValue=-3.4028235e+38 $gmt -version a.nc.zip

welcome[bot] commented 2 weeks ago

👋 Thanks for opening your first issue here! Please make sure you filled out the template with as much detail as possible. We appreciate that you took the time to contribute!

Please make sure you read our Contributing Guide and abide by our Code of Conduct.

joa-quim commented 2 weeks ago

It doesn't report them because they are not in the grid's metadata (as they should). Run

gmt grdinfo a.nc -L

and you'll see them.