Opened 4 years ago

Closed 4 years ago

Last modified 4 years ago

#5118 closed defect (fixed)

NetCDF driver doesn't recognize grids as uniform.

Reported by: Kyle Shannon Owned by: etourigny
Priority: normal Milestone: 1.10.1
Component: default Version: unspecified
Severity: normal Keywords: netcdf
Cc: etourigny

Description

In GDAL 1.10.x the checks for uniform grids use 5 decimal places. Prior to 1.10.x, it used 2 decimal places. While increased precision is probably the correct way to handle the situation, it doesn't work on grids it used to work on. The attached patch just allows a config option to revert back to the old behavior using NETCDF_SPACING_PRECISION. Etienne, let me know if you think this is worth adding. I would just like to be able to treat files that worked in 1.9 the same in 1.10, even though it's not totally correct.

Attachments (4)

netcdf_space.diff (3.3 KB) - added by Kyle Shannon 4 years ago.
20130614T2100.nc (109.9 KB) - added by Kyle Shannon 4 years ago.
5118_2.txt (5.7 KB) - added by etourigny 4 years ago.
5118_3.txt (12.1 KB) - added by etourigny 4 years ago.

Download all attachments as: .zip

Change History (23)

Changed 4 years ago by Kyle Shannon

Attachment: netcdf_space.diff added

comment:1 Changed 4 years ago by etourigny

Owner: changed from warmerdam to etourigny

Kyle - I was afraid this change might lead to some issues. Can you please upload the file(s) which cause the trouble? or maybe just the lat/lon grids.

comment:2 Changed 4 years ago by etourigny

actually, before it used 4 decimal places.

comment:3 Changed 4 years ago by etourigny

I'm not sure a new config option is the way to go, ideally the code should handle all situations. I tried using CplIsEqual?() instead of this method but it still failed...

comment:4 Changed 4 years ago by etourigny

see related bug #5114 which prompted this change

Changed 4 years ago by Kyle Shannon

Attachment: 20130614T2100.nc added

comment:5 Changed 4 years ago by Kyle Shannon

Etienne, the grid is more non-uniform than I expected. I originally just looked at ncdump. Inspecting the values in a debugger showed that the grids are farther off than what I expected.

ncdump:

kyle@kyle-workstation:build$ ncdump -v x ../trunk/data/NCEP-RAP-13km-SURFACE-big_butte.tif/20130614T2100/20130614T2100.nc 
netcdf \20130614T2100 {
dimensions:
	time = 17 ;
	height_above_ground = 1 ;
	y = 22 ;
	x = 18 ;
	height_above_ground1 = 1 ;
variables:
...
data:

 x = -1652.575, -1639.03, -1625.485, -1611.94, -1598.395, -1584.85, 
    -1571.305, -1557.76, -1544.215, -1530.67, -1517.125, -1503.58, -1490.035, 
    -1476.49, -1462.945, -1449.4, -1435.855, -1422.31 ;

gdb:

(gdb) p pdfXCoord[0]@18
$1 = {-1652.5753173828125, -1639.0302734375, -1625.4853515625, -1611.9403076171875, -1598.3953857421875, -1584.850341796875, -1571.3052978515625, -1557.7603759765625, -1544.21533203125, 
  -1530.6702880859375, -1517.1253662109375, -1503.580322265625, -1490.0352783203125, -1476.4903564453125, -1462.9453125, -1449.4002685546875, -1435.8553466796875, -1422.310302734375}

I believe that the people producing the grids have some issues.

kss

comment:6 Changed 4 years ago by etourigny

Kyle,

I have a feeling these issues are due to errors in rounding, when the lon/lat or x/y variables are defined as float. Because if you look at the values given by ncdump (output as float), the spacing is uniform, whereas in gdal it is not (output as double)- because of rounding errors in float->double conversion. I think the code should use floats when the dimension variables are floats, and check for diff<float precision.

Also, the 0.1 check if for gaussian grids, which should only be used when the variable is geographic (not projected), what do you think?

comment:7 in reply to:  6 Changed 4 years ago by Kyle Shannon

Replying to etourigny:

Kyle,

I have a feeling these issues are due to errors in rounding, when the lon/lat or x/y variables are defined as float. Because if you look at the values given by ncdump (output as float), the spacing is uniform, whereas in gdal it is not (output as double)- because of rounding errors in float->double conversion. I think the code should use floats when the dimension variables are floats, and check for diff<float precision.

It could be, but it seems like the decimal place of the precision isn't that far out.

Also, the 0.1 check if for gaussian grids, which should only be used when the variable is geographic (not projected), what do you think?

This is probably true too, but I don't have much of an opinion. It seems correct.

kss

comment:8 Changed 4 years ago by etourigny

in this case, it seems that using single instead of double doesn't work:

there is a difference of 0.000122 between the spacing at begin and at middle, which is larger that single-point precision.

THere needs to have a comparison that scales with the amplitude of the data.

Output from some experimental code:

GDAL_netCDF: dfSpacingBegin: 13.545044 dfSpacingMiddle: 13.544922 dfSpacingLast: 13.545044
GDAL_netCDF: dfSpacingBegin: 13.545 dfSpacingMiddle: 13.5449 dfSpacingLast: 13.545
GDAL_netCDF: fSpacingBegin: 13.545044 fSpacingMiddle: 13.544922 fSpacingLast: 13.545044
GDAL_netCDF: fSpacingBegin: 13.545 fSpacingMiddle: 13.5449 fSpacingLast: 13.545
GDAL_netCDF: xdim: 18 nSpacingBegin: 13545 nSpacingMiddle: -13545 nSpacingLast: -13545
GDAL_netCDF: xcoords: -1652.575317 -1639.030273 -1530.670288 -1517.125366 -1435.855347 -1422.310303
GDAL_netCDF: xdiff: -13.545044 -13.544922 -13.545044 / 0.000122 0.000000
GDAL_netCDF: xcoords float: -1652.58 -1639.03 -1530.67 -1517.13 -1435.86 -1422.31
GDAL_netCDF: xdiff float: -13.545044 -13.544922 -13.545044 / 0.000122 0.000000

comment:9 Changed 4 years ago by Kyle Shannon

Etienne, I dug into this a little more, and I my fix isn't even remotely close. I tried to warp it and it fails. I am not familiar with the geolocation arrays interface, but I would think that it needs x and y coordinates. In this case, should geolocation arrays be used?

kyle@kyle-workstation:20130618T2100$ CPL_DEBUG=GDAL_netCDF gdalwarp -t_srs epsg:32612 NETCDF:"20130618T2100.nc":v-component_of_wind_height_above_ground test.tif -tr 12000 12000
GDAL_netCDF: 
=====
Open(), filename=[test.tif]
GDAL_netCDF: 
=====
Open(), filename=[test.tif]
GDAL_netCDF: 
=====
Open(), filename=[NETCDF:20130618T2100.nc:v-component_of_wind_height_above_ground]
GDAL_netCDF: calling nc_open( 20130618T2100.nc )
GDAL_netCDF: got cdfid=65536

GDAL_netCDF: driver detected file type=1, libnetcdf detected type=1
GDAL_netCDF: dim_count = 5
GDAL_netCDF: var_count = 10
GDAL_netCDF: variable #1 [time] was ignored
GDAL_netCDF: variable #2 [height_above_ground] was ignored
GDAL_netCDF: variable #3 [y] was ignored
GDAL_netCDF: variable #4 [x] was ignored
GDAL_netCDF: variable #7 [height_above_ground1] was ignored
GDAL_netCDF: 
=====
SetProjectionFromVar( 6 )

GDAL_netCDF: got grid_mapping LambertConformal_Projection
GDAL_netCDF: bIsGdalFile=0 bIsGdalCfFile=0 bBottomUp=1
GDAL_netCDF: got spheroid from CF: (6371229.000000 , -1.000000)
GDAL_netCDF: set bBottomUp = 1 from Y axis
GDAL_netCDF: xdim: 18 nSpacingBegin: 135449 nSpacingMiddle: -135450 nSpacingLast: -135450
GDAL_netCDF: Longitude is not equally spaced.
GDAL_netCDF: ydim: 23 nSpacingBegin: 135452 nSpacingMiddle: -135449 nSpacingLast: -135449
Warning 1: Latitude grid not spaced evenly.
Seting projection for grid spacing is within 0.1 degrees threshold.

GDAL_netCDF: Latitude grid not spaced evenly, but within 0.1 degree threshold (probably a Gaussian grid).
Saving original latitude values in Y_VALUES geolocation metadata
GDAL_netCDF: units=km
GDAL_netCDF: setting WKT from CF
GDAL_netCDF: SetProjection, WKT = PROJCS["unnamed",GEOGCS["unknown",DATUM["unknown",SPHEROID["Sphere",6371229,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_1SP"],PARAMETER["latitude_of_origin",25],PARAMETER["central_meridian",265],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],PARAMETER["standard_parallel_1",25],UNIT["kilometre",1000,AUTHORITY["EPSG","9036"]]]
GDAL_netCDF: got SRS but no geotransform from CF!
GDAL_netCDF: coordinates attribute [time height_above_ground1 y x ] is unsupported
GDAL_netCDF: bGotGeogCS=1 bGotCfSRS=1 bGotCfGT=0 bGotGdalSRS=0 bGotGdalGT=0
GDAL_netCDF: did not get geotransform from CF nor GDAL!
GDAL_netCDF: netcdf type=5 gdal type=6 signedByte=1
GDAL_netCDF: netcdf type=5 gdal type=6 signedByte=1
GDAL_netCDF: netcdf type=5 gdal type=6 signedByte=1
GDAL: GDALOpen(NETCDF:20130618T2100.nc:v-component_of_wind_height_above_ground, this=0x20220f0) succeeds as netCDF.
ERROR 1: Missing some geolocation fields in GDALCreateGeoLocTransformer()
GDAL: GDALClose(NETCDF:20130618T2100.nc:v-component_of_wind_height_above_ground, this=0x20220f0)
GDAL: In GDALDestroy - unloading GDAL shared library.

comment:10 Changed 4 years ago by etourigny

1D geolocation arrays are supported, when e.g. one dimension does not have fixed values. I guess something is missing, I will try to look into this. I implemented this some time ago and only partially tested...

What do you suggest we do here? Should the x/y spacing be treated as uniform or no?

comment:11 Changed 4 years ago by etourigny

Looking back into the code... the Y_VALUES geoloction metadata is a kind of hack I added to be able to preserve the original y values in the case of gaussian grids (within 0.1 deg threshold). A side-effect of this is that gdalwarp is trying to use the (fake) geolocation info (because there is a geoloation metadata domain), but it fails to find the necessary info like X_BAND and Y_BAND, which do not exist.

It would be pretty involved to fix this... perhaps the Y_VALUES metadata should be placed elsewhere than in the Geolocation domain.

Changed 4 years ago by etourigny

Attachment: 5118_2.txt added

comment:12 Changed 4 years ago by Kyle Shannon

Etienne, There are a few ways to look at it. The grids aren't uniform, so they maybe shouldn't be treated as so. The check for the gaussian stuff is in degrees. The file I am looking at is projected, and the values are in kilometers when they are checked and the difference is on the order of meters, which is probably reasonable, but maybe a config option would be appropriate. I have no strong feeling either way, other than I can't read/warp these files anymore.

If the Y_VALUES aren't actual Geolocation arrays, they probably shouldn't be stored in the Geolocation arrays.

comment:13 Changed 4 years ago by etourigny

Kyle,

if you look at the data with ncdump, the data IS uniform. It seems to me that the difference is because of a conversion from float to double...

Can you test that it works with the supplied patch 5118_2.txt? If so I will commit because I see this as a regression.

It could be valid to add the "gaussian" check only for geographic (lat/lon) grids, and not for projected grids.

Also, I will put the Y_VALUES in another domain, say geolocation2.

Changed 4 years ago by etourigny

Attachment: 5118_3.txt added

comment:14 Changed 4 years ago by etourigny

Please review patch 5118_3.txt which resolves all issues, hopefully.

comment:15 Changed 4 years ago by Kyle Shannon

Etienne, The patch fixes everything for me. I'm okay with it if you are. The metadata domain fix is good too, the warper won't freak out on it. Thanks.

comment:16 Changed 4 years ago by Kyle Shannon

Thanks a lot Etienne,

kss

comment:17 Changed 4 years ago by etourigny

added autotest to trunk in r26107 and fixes in r26108 and r26109

comment:18 Changed 4 years ago by etourigny

Resolution: fixed
Status: newclosed

fixed in 1.10 (r26110)

comment:19 in reply to:  16 Changed 4 years ago by etourigny

Replying to kyle:

Thanks a lot Etienne,

kss

thanks for your help

Note: See TracTickets for help on using tickets.