#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)
Change History (23)
by , 11 years ago
Attachment: | netcdf_space.diff added |
---|
comment:1 by , 11 years ago
Owner: | changed from | to
---|
comment:3 by , 11 years ago
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...
by , 11 years ago
Attachment: | 20130614T2100.nc added |
---|
comment:5 by , 11 years ago
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
follow-up: 7 comment:6 by , 11 years ago
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 by , 11 years ago
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 by , 11 years ago
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 by , 11 years ago
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 by , 11 years ago
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 by , 11 years ago
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.
by , 11 years ago
Attachment: | 5118_2.txt added |
---|
comment:12 by , 11 years ago
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 by , 11 years ago
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.
by , 11 years ago
Attachment: | 5118_3.txt added |
---|
comment:15 by , 11 years ago
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 - 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.