Opened 17 years ago
Last modified 17 years ago
#1506 closed defect
gdalinfo on a netCDF file returns incorrect coordinates — at Version 13
Reported by: | Owned by: | ||
---|---|---|---|
Priority: | normal | Milestone: | 1.5.0 |
Component: | GDAL_Raster | Version: | 1.4.0 |
Severity: | normal | Keywords: | NetCDF |
Cc: | dem@…, warmerdam, Mateusz Łoskot |
Description (last modified by )
The attached netCDF file was generated using NCO (ncra -- see "history"), which has the dimensions: lon = 128 ; lat = 64 ; time = UNLIMITED ; (1 currently). This is effectively a raster showing average global temperature (the time was averaged between 2000-2100).
However, using "gdalinfo" 1.4.0 on Cygwin WinXPSP2:
$ gdalinfo out.nc Driver: netCDF/Network Common Data Format Size is 512, 512 Coordinate System is `' Metadata: NC_GLOBAL#title=CCCma model output prepared for IPCC Fourth Assessment 720 ppm stabilization experiment (SRES A1B) NC_GLOBAL#institution=CCCma (Canadian Centre for Climate Modelling and Analysis, Victoria, BC, Canada) NC_GLOBAL#source=CGCM3.1 (2004): atmosphere: AGCM3 (GCM13d, T63L31); ocean: CCCMA (OGCM3.1,256x192L29) NC_GLOBAL#contact=Greg Flato (Greg.Flato@ec.gc.ca) NC_GLOBAL#project_id=IPCC Fourth Assessment NC_GLOBAL#table_id=Table A1 (17 November 2004) NC_GLOBAL#experiment_id=720 ppm stabilization experiment (SRES A1B) NC_GLOBAL#realization=1 NC_GLOBAL#cmor_version=9.600000e-01 NC_GLOBAL#Conventions=CF-1.0 NC_GLOBAL#history=Wed Feb 28 16:27:16 2007: ncra -v ts ts_a1_sresa1b_1_cgcm3.1_t63_2001_2100.nc out.nc At 00:07:46 on 06/21/2005, CMOR rewrote data to comply with CF standards and IPCC Fourth Assessment requirements NC_GLOBAL#comment=This model run continues from the end of the 20th century simulation with GHG and aerosol loadings for the IPCC SRES A1B scenario as tabulated in the IPCC Third Assessment Report, Appendix II. The CO2 concentrations are from the Bern-CC model (Contribution of Working Group I to the Third Assesment Report of IPCC, p808) and the aerosol loadings are from O. Boucher, Laboratoire d'Optique Atmospherique, France. For years 2101-2300, all GHG concentrations and the aerosol loading are held constant at the values obtained by extrapolation to year 2101. NC_GLOBAL#nco_openmp_thread_number=1 Subdatasets: SUBDATASET_1_NAME=NETCDF:"out.nc":ts SUBDATASET_1_DESC=[1x64x128] surface_temperature (32-bit floating-point) Corner Coordinates: Upper Left ( 0.0, 0.0) Lower Left ( 0.0, 512.0) Upper Right ( 512.0, 0.0) Lower Right ( 512.0, 512.0) Center ( 256.0, 256.0)
I've experienced very similar behaviour with other netCDF files, which have very different dimensions, but still report 512x512 using gdalinfo.
Change History (15)
by , 17 years ago
comment:1 by , 17 years ago
Cc: | added |
---|---|
Description: | modified (diff) |
Milestone: | → 1.4.1 |
Owner: | changed from | to
Julien Demaria reports the following via the mailing list, and this is presumably related. I have reassigned this bug report to Denis who is responsible for this driver. It would be really nice to have this resolved for GDAL 1.4.1.
Hi,
I have found a bug in the computing of the GeoTransform in the function netCDFDataset::SetProjection of the netcdfdataset.cpp file (last version): the resolution was wrong and the origin was the center and not the upper-left corner:
poDS->adfGeoTransform[0] = pdfXCoord[0]; poDS->adfGeoTransform[3] = pdfYCoord[0]; poDS->adfGeoTransform[2] = 0;
poDS->adfGeoTransform[4] = 0; poDS->adfGeoTransform[1] = (( pdfXCoord[xdim-1] -
pdfXCoord[0] ) /
poDS->nRasterXSize);
poDS->adfGeoTransform[5] = (( pdfYCoord[ydim-1] -
pdfYCoord[0] ) /
poDS->nRasterYSize);
should be corrected with:
poDS->adfGeoTransform[0] = pdfXCoord[0]; poDS->adfGeoTransform[3] = pdfYCoord[0]; poDS->adfGeoTransform[2] = 0;
poDS->adfGeoTransform[4] = 0; poDS->adfGeoTransform[1] = (( pdfXCoord[xdim-1] -
pdfXCoord[0] ) /
( poDS->nRasterXSize - 1 ));
poDS->adfGeoTransform[5] = (( pdfYCoord[ydim-1] -
pdfYCoord[0] ) /
( poDS->nRasterYSize - 1 ));
- (poDS->adfGeoTransform[1] / 2);
- (poDS->adfGeoTransform[5] / 2);
Maybe there is also the same problem just after this code when trying to retrieve the geolocation from the attributes like "Northernmost_Northing", depending if these attributes define border or center pixel:
adfGeoTransform[0] = dfWE; adfGeoTransform[1] = (dfEE - dfWE) / poDS->GetRasterXSize(); adfGeoTransform[2] = 0.0; adfGeoTransform[3] = dfNN; adfGeoTransform[4] = 0.0; adfGeoTransform[5] = (dfSN - dfNN) / poDS->GetRasterYSize();
Regards,
Julien
follow-up: 5 comment:3 by , 17 years ago
Component: | default → GDAL_Raster |
---|---|
Resolution: | → fixed |
Status: | new → closed |
Hi Julien,
Thank you for your input. There was effectively a computation error in Set Projection. I have committed your change to compute the pixel center.
If you see other errors like this one, do not hesitate to communicate them to me and I will make the appropriate changes.
BTW, in the attached file (out.nc) the lat/lon variable do not contain data equally space. (Is this correct?)
The netcdf driver does not handle non-gridded data....
Best Regards, Denis
comment:4 by , 17 years ago
I have applied this change to the 1.4 branch as well, so it will make it into 1.4.1 (r11169).
comment:5 by , 17 years ago
Hi Denis, Great work, I'll have to test your commits later today.
BTW, in the attached file (out.nc) the lat/lon variable do not contain data equally space. (Is this correct?)
Regarding the "out.nc" file: this particular file has equal spaced longitudes, and non-equal latitudes, which is typical of Gaussian grids (see http://en.wikipedia.org/wiki/Gaussian_grid for more info). For the purposes of GDAL (and due to practical limitations in the architecture), this can be generalized by single 'dx' and 'dy' values.
comment:6 by , 17 years ago
Resolution: | fixed |
---|---|
Status: | closed → reopened |
I have compiled GDAL (r11169) on Debian Sid (linux 2.6.18-4-686), however gdalinfo reports the exact (and incorrect) coordinate information for out.nc (as well as other nc files I use). Has the error also been addressed in the gdalinfo command?
comment:7 by , 17 years ago
(sorry, I meant to say "exact same" incorrect coordinate information as before; e.g., 0.0, 256.0, 512, etc.)
comment:8 by , 17 years ago
I think this bug contains 2 different topics.
First, I have commit changes written by Julien Demaria to set the projection to the center of a pixel instead of its border.
Second, I think what mwtoews is saying is that no projection is set for the attached file out.nc. Since this file is not equally spaced in longitude and latitude this is normal. GDAL netCDF driver is limited to equally spaced lat/lon.
I have looked at the Gaussian_grid web page and found out that this projection model is not supported by proj4 either.
This grid cannot be handled by GDAL GeoTransform.
Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2) Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)
Therefore I do not know how to set the projection for this type of gridded data.
For the attached file 'out.nc' the following results are correct and expected by gdalinfo.
gdalinfo -nomd ./out.nc Driver: netCDF/Network Common Data Format Size is 512, 512 Coordinate System is `' Subdatasets: SUBDATASET_1_NAME=NETCDF:"./out.nc":ts SUBDATASET_1_DESC=[1x64x128] surface_temperature (32-bit floating-point) Corner Coordinates: Upper Left ( 0.0, 0.0) Lower Left ( 0.0, 512.0) Upper Right ( 512.0, 0.0) Lower Right ( 512.0, 512.0) Center ( 256.0, 256.0) gdalinfo -nomd NETCDF:"./out.nc":ts Driver: netCDF/Network Common Data Format Size is 128, 64 Coordinate System is `' Corner Coordinates: Upper Left ( 0.0, 0.0) Lower Left ( 0.0, 64.0) Upper Right ( 128.0, 0.0) Lower Right ( 128.0, 64.0) Center ( 64.0, 32.0) Band 1 Block=128x1 Type=Float32, ColorInterp=Undefined NoData Value=9.96920996838687e+36
Denis
comment:9 by , 17 years ago
Denis, I'll admit that the 'out.nc' file I uploaded was an odd example (I was trying to illustrate a simple example which displays the strange power of 2 coordinate information). I don't think it is possible to define a projection for a Gaussian grid (it depends on the Gaussian quadrature for the wave number used for spectral finite differences).
For the purposes of the GDAL driver, these types of netCDF grids (and they are common) can be assumed to have equally spaced latitudes. For this particular file, the mean distances between the grid points in the latitude dimension is 2.789327º, so the error in that assumption is between -0.029° and 0.002°, which is very minor (considering the scale). The outer edge of latitude would then be +/-89.25846º, which is within +/-90º.
Thanks for alerting me to the correct usage for 'gdalinfo'; I was unaware that I needed to specify the variable within the command. (This does still seem a bit odd since netCDF files only have a single set of latitude/longitude dimensions, however I see the complexity involved since each variable can be defined using different dimension variables. I would rather see "undefined" or an equivalent, rather than false "Size is 512, 512" information.)
follow-up: 11 comment:10 by , 17 years ago
Michael,
I could set a threshold for latitude for this kind of lat/lon data. People who have a different projection (UTM,LCC) will not get into this kind of problem. (Lat/Lon data usually have a similar scale as out.nc.)
Frank, do you think I should display a warning when the data are not equally spaced, but are withing an acceptable threshold for "Set projection" method? I will have to document this behavior.
Michael, as for specifying the variable withing the 'gdalinfo' command, normaly when the file contains only a single variable you should not have to specify that variable. This is something I never finished getting around. (Now that you got caught by this problem, it will push me to update my code...)
As for specifying "undefined" instead of a false size, this is an issue that should be address to Frank Warmerdam as well.
comment:11 by , 17 years ago
Replying to dnadeau:
Frank, do you think I should display a warning when the data are not equally spaced, but are withing an acceptable threshold for "Set projection" method? I will have to document this behavior.
Denis,
This approach is fine with me, though I consider this a 1.5.0 issue.
Michael, as for specifying the variable withing the 'gdalinfo' command, normaly when the file contains only a single variable you should not have to specify that variable. This is something I never finished getting around. (Now that you got caught by this problem, it will push me to update my code...)
I'd dearly like this aspect finished up!
As for specifying "undefined" instead of a false size, this is an issue that should be address to Frank Warmerdam as well.
It is standard to use a default 512x512 with no bands to represent a strictly container oriented holder of subdatasets. This is not unique to this driver.
BTW, if the core bug is fixed, please re-milestone this ticket to 1.5.0
comment:12 by , 17 years ago
The issues described by Julien and myself in #1505 in regards to the centre/edge issue are resolved from r11169 and r11175; so the core bug is resolved.
The issue at hand (possibly for a later milestone) is now in regards to non-equally-spaced latitude grid-points, as used in Gaussian grids. I'm not involved in any way in regards to creating these data sets, so I'm by no means an expert; however, I am partially familiar with them, since I use their data. A good reference and resource on GCM grid types is http://www.ncl.ucar.edu/Document/Functions/sphpk_grids.shtml.
In my opinion, I think a warning is appropriate for all grids with non-equally-spaced latitudes, but proceed with a "best guess" latitude spacing from the average. There are tools and resources for scientists to properly regrid between fixed and Gaussian grids, for example, the tools in NCL:
However, if you decide to proceed with determining an acceptable threshold for the warning, I have provided a plot (attached to this ticket). This plot was created form the latitude dimensions of two GCMs from CCCma, and from the NCEP/NCAR Reanalysis Project. I calculated the normalized error here based on (using R syntax with a file from NCEP/NCAR):
nc <- open.ncdf("air2m.mon.mean.nc") lat <- get.var.ncdf(nc,"lat") dlat <- abs(diff(lat)) err <- (dlat/mean(dlat)-1)*100
From my plot, the NCEP/NCAR normalized error is the largest at the two ends at +/-0.799749978. Therefore, an acceptable threshold should be +/-0.8 % of the normalized latitude differences. However, as my opinion above states, I think a warning should be thrown regardless, as there are good tools external to GDAL to regrid this kind of data.
comment:13 by , 17 years ago
Description: | modified (diff) |
---|
Example netCDF file