gdal_translate misinterprets netCDF dimensions and cellsize

[mwtoews@…]

This issue was initially posted on the gdal-dev mailing list using GDAL 1.3.2 (see ​http://lists.maptools.org/pipermail/gdal-dev/2006-October/010393.html for original message and netCDF file attachment). I am experiencing the same behaviour in GDAL 1.4.0 under Cygwin WinXPSP2.

My netCDF file was generated using NCO (ncra), and now I need to convert the NetCDF file into a different format for GIS programs:

$ gdal_translate -of AAIGrid -sds test.nc test.asc
$ mv test.asc1 test.asc

No errors were generated during the conversion.

The header of the AAIGrid output is:

ncols        86
nrows        225
xllcorner    747700.000000000000
yllcorner    5441461.627906977199
dx           98.837209302326
dy           99.555555555556
NODATA_value  -9999

However, these headers should read:

ncols         86
nrows         225
xllcorner     747700
yllcorner     5441200
cellsize      100
NODATA_value  -9999

The dx/dy in the output of the file was incorrectly calculated as:

dx=cellsize*(ncol-1)/ncol
dy=cellsize*(nrow-1)/nrow

theses formulas are of course nonsense, but are apparently how they were calculated. This is also why dx != dy (because ncols != nrows), and why they are not correct, since they are multiplied by the ratio of base0/base1 representations of their dimension lengths.

My netCDF file (attached on original forum post -- link at top of bug report) should have a cellsize of exactly 100, since my data are spaced every 100 m in the x- and y-dims. This netCDF file was produced using a custom python script, which uses Scientific.IO.NetCDF for python 2.3. The file was later averaged, using 'ncra', which is part of NCO (command viewed in history' global attribute). My netCDF file is not perfect; I've recently discovered that my 'ydim' is 5463700, 5463600, ..., 5441300; or it is in reverse order (dy = -100). For netCDF files, this should not matter (although I'll admit that it is a bit sloppy); I've been able to use the file as normal with 'ncview' and NCO.

(As a caution, I have also seen some netCDF files which order their latitudes from North to South, which is also in reverse order; I'm not sure if this is tested, and corrected in GDAL, but I would highly advise it for netCDF files).

[warmerdam]

Assigning to Denis under the assumption this is a netcdf problem, and likely related to ticket #1506. Also marked as 1.4.1 milestone if possible.

[dnadeau]

I have downloaded the file from the forum and I got different results from you.

ncols        86
nrows        225
xllcorner    747650.000000000000
yllcorner    5441250.000000000000
cellsize     100.000000000000
NODATA_value  -9999

The corners differ from 50 units for GDAL netCDF driver now derives the center of the pixel.

As well, I could not find the following "nonsense" formula in the netcdfdataset.cpp file. Can you be more specific, I think they should be deal with right away? dx=cellsize*(ncol-1)/ncol dy=cellsize*(nrow-1)/nrow

What version of GDAL do you have?

If you can provide me with netCDF files which order their latitudes from North to South, I will be very happy to test them out.

Regards, Denis

Denis

[Mike Taves]

The "nonsense" formula was just a back-calculation of the results from gdal_translate; I have not looked at netcdfdataset.cpp for this. The issue was in versions 1.3.2 and 1.4.0, and can be reproduced in a FWTools 1.2.2 Shell, which has GDAL 14.0.0 released 2007/02/22. The dx and dy cellsize issue appears to be resolved from your most recent revisions related to #1506.

However, there is a new (but possibly related) issue: I have just downloaded GDAL through SVN (r11169) and compiled GDAL in Cygwin, and I see (using the same file/conversion):

ncols        86
nrows        225
xllcorner    37385000.000000000000
yllcorner    -273207500.000000000000
cellsize     100.000000000000
NODATA_value  -9999

Interestingly, the [xy]llcorner reported here are approximately 50 and -50 times what you have reported; or in other words, half of the dimensioned cellsize.

This version of GDAL is compiled with netcdf library version 3.6.1-beta3 of Feb 17 2006 (the same netcdf library I used for my initial report). Also, the version from gdal_translate.exe is reported to be "GDAL 1.4.0.0, released 2007/01/02", even though I just downloaded/compiled it (shouldn't this say differently?).

I have also just compiled netCDF 3.6.2 and GDAL (r11169) on Debian Sid, and can confirm your correct results (also the version reports GDAL 1.5dev 2007/04/01).

The test-0001.nc file is already ordered from North to South (using Northings, not Latitudes) -- do you mean positive increasing order, from South to North? I plan to have this resolved with my setup in a few days, and will attach a similar file here.

And thanks for the correction and clarification for the corner positions.

[dnadeau]

Michael,

Try to update netcdfdataset.cpp (r11175). When I cut/paste my code, a substraction sign was converted into a multiplication when computing xllcorner. (wiki formatting...)

[Mike Taves]

That must have fixed it! Both my Debian sid and cygwin versions of GDAL (r11175) now show the same and correct coordinate information from the translation! Thanks for your fix!