NetCDF driver doesn't recognize grids as uniform.

[Kyle Shannon]

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.

[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.

[etourigny]

actually, before it used 4 decimal places.

[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...

[etourigny]

see related bug #5114 which prompted this change

[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

[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?

[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

[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

[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.

[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?

[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.

[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.

[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.

[etourigny]

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

[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.

[Kyle Shannon]

Thanks a lot Etienne,

kss

[etourigny]

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

[etourigny]

fixed in 1.10 (r26110)

[etourigny]

Replying to kyle:

Thanks a lot Etienne,

kss

thanks for your help