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

[warmerdam]

I imagine ticket #1505 is the same problem.

[dnadeau]

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

[warmerdam]

I have applied this change to the 1.4 branch as well, so it will make it into 1.4.1 (r11169).

[Mike Taves]

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.

[Mike Taves]

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?

[Mike Taves]

(sorry, I meant to say "exact same" incorrect coordinate information as before; e.g., 0.0, 256.0, 512, etc.)

[dnadeau]

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

[Mike Taves]

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

[dnadeau]

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.

[warmerdam]

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

[Mike Taves]

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:

• ​g2fsh - from a Gaussian grid to a fixed grid
• ​f2gsh - from a fixed grid to a Gaussian grid

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.

[Mike Taves]

Errors in Gaussian grid latitude spacings

[warmerdam]

Remaining issue deferred to 1.5.0. It likely ought to be split out into a new focused ticket, and this ticket closed.

[dnadeau]

Michael,

I have changed the netcdfdataset to reflect our comments.

Could you please validate the following information from you file? I will then check it in.

Thanks in advance.

Regards, Denis

gdalinfo NETCDF:"/disk3/dnadeau/download/out.nc":ts Warning 1: Latitude grid not spaced evenly. Seting projection for grid spacing is within 0.1 degrees threshold.

Driver: netCDF/Network Common Data Format Size is 128, 64 Coordinate System is `' Origin = (-1.406250000000000,-89.258462312871046) Pixel Size = (2.812500000000000,2.789326947277220) 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@…) 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 ts#standard_name=surface_temperature ts#long_name=Surface Skin Temperature ts#units=K ts#cell_methods=time: mean (interval: 15 minutes) ts#original_name=GT + 273.16 ts#history= At 00:07:46 on 06/21/2005: CMOR altered the data in the following ways: converted from type "double" to type "real"; lon#standard_name=longitude lon#long_name=longitude lon#units=degrees_east lon#axis=X lon#bounds=lon_bnds lat#standard_name=latitude lat#long_name=latitude lat#units=degrees_north lat#axis=Y lat#bounds=lat_bnds time#standard_name=time time#long_name=time time#units=days since 1850-1-1 time#axis=T time#calendar=365_day time#bounds=time_bnds

Corner Coordinates: Upper Left ( -1.4062500, -89.2584623) Lower Left ( -1.4062500, 89.2584623) Upper Right ( 358.594, -89.258) Lower Right ( 358.594, 89.258) Center ( 178.5937500, 0.0000000) Band 1 Block=128x1 Type=Float32, ColorInterp=Undefined

NoData Value=9.96920996838687e+36 Metadata:

NETCDF_VARNAME=ts NETCDF_DIMENSION_time=73364.2 NETCDF_time_units=days since 1850-1-1

[Mike Taves]

Denis,

Those numbers look good (given the info). I checked the calculations in R:

library(ncdf)
nc <- open.ncdf("out.nc")
lat <- get.var.ncdf(nc,"lat")
dlat <- diff(lat)
dy <- mean(dlat)
lon <- get.var.ncdf(nc,"lon")
dlon <- diff(lon)
dx <- mean(dlon)

> dx # x-dim pixel size
[1] 2.8125
> dy # mean x-dim pixel size
[1] 2.789327
> lat[1]-dy/2 # bottom edge
[1] -89.25846
> lat[length(lat)]+dy/2 # top edge
[1] 89.25846
> lon[1]-dx/2 # left edge
[1] -1.40625
> lon[length(lon)]+dx/2 # right edge
[1] 358.5938
> (lon[length(lon)]-lon[1])/2 # center in x-dim
[1] 178.5938
> (lat[1]+lat[length(lat)])/2 # center in y-dim
[1] 0

And the warning looks appropriate. You can close this ticket after your commits, and I'll refile a ticket for the remaining issue about showing coordinate information for a netCDF file with only one variable, as you mentioned earlier. Thanks!

[dnadeau]

Threshold has been comitted (r11237).