gdalinfo on a netCDF file returns incorrect coordinates

[mwtoews@…]

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.

[mwtoews@…]

Example netCDF file

[warmerdam]

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[0] = pdfXCoord[0]

• (poDS->adfGeoTransform[1] / 2);

poDS->adfGeoTransform[3] = pdfYCoord[0]

• (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