NetCDF driver writes incorrect GeoTranform in case of WRITE_BOTTOMUP=YES

[korosov]

In both cases when I use WRITE_BOTTOMUP=YES or WRITE_BOTTOMUP=NO the last term of GeoTransform in metadata is negative. However, in the case of WRITE_BOTTOMUP=YES the Y coordinate is increasing, therefore the increment in the GeoTransform should be positive. Attached are two NC files generated with WRITE_BOTTOMUP=YES and WRITE_BOTTOMUP=NO using the following commands:

gdalwarp -overwrite -of NetCDF -t_srs '+proj=latlong +no_defs' -te 27 70 31 72 -ts 50 50 -co 'WRITE_BOTTOMUP=NO' gcpb1.tif gcps_pro_no.nc

gdalwarp -overwrite -of NetCDF -t_srs '+proj=latlong +no_defs' -te 27 70 31 72 -ts 50 50 -co 'WRITE_BOTTOMUP=YES' gcpb1.tif gcps_pro_yes.nc

ncdump -h gives:

netcdf gcps_pro_no { dimensions:

lon = 50 ; lat = 50 ;

variables:

byte Band1(lat, lon) ;

Band1:long_name = "GDAL Band Number 1" ; Band1:_Unsigned = "true" ; Band1:valid_range = 0s, 255s ; Band1:_FillValue = 0b ; Band1:grid_mapping = "crs" ;

char crs ;

crs:grid_mapping_name = "latitude_longitude" ; crs:longitude_of_prime_meridian = 0. ; crs:semi_major_axis = 6378137. ; crs:inverse_flattening = 298.257223563 ; crs:spatial_ref = "GEOGCS[\"unnamed ellipse\",DATUM[\"unknown\",SPHEROID[\"unnamed\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]]" ; crs:GeoTransform = "27 0.08 0 72 0 -0.04 " ;

double lat(lat) ;

lat:standard_name = "latitude" ; lat:long_name = "latitude" ; lat:units = "degrees_north" ;

double lon(lon) ;

lon:standard_name = "longitude" ; lon:long_name = "longitude" ; lon:units = "degrees_east" ;

global attributes:

:Conventions = "CF-1.5" ; :GDAL = "GDAL 2.0dev, released 2011/12/29" ; :history = "Tue Nov 06 10:07:11 2012: GDAL Create( gcps_pro_no.nc, ... )" ;

data:

lat = 71.98, 71.94, 71.9, 71.86, 71.82, 71.78, 71.74, 71.7, 71.66, 71.62,

71.58, 71.54, 71.5, 71.46, 71.42, 71.38, 71.34, 71.3, 71.26, 71.22, 71.18, 71.14, 71.1, 71.06, 71.02, 70.98, 70.94, 70.9, 70.86, 70.82, 70.78, 70.74, 70.7, 70.66, 70.62, 70.58, 70.54, 70.5, 70.46, 70.42, 70.38, 70.34, 70.3, 70.26, 70.22, 70.18, 70.14, 70.1, 70.06, 70.02 ;

lon = 27.04, 27.12, 27.2, 27.28, 27.36, 27.44, 27.52, 27.6, 27.68, 27.76,

27.84, 27.92, 28, 28.08, 28.16, 28.24, 28.32, 28.4, 28.48, 28.56, 28.64, 28.72, 28.8, 28.88, 28.96, 29.04, 29.12, 29.2, 29.28, 29.36, 29.44, 29.52, 29.6, 29.68, 29.76, 29.84, 29.92, 30, 30.08, 30.16, 30.24, 30.32, 30.4, 30.48, 30.56, 30.64, 30.72, 30.8, 30.88, 30.96 ;

}

netcdf gcps_pro_yes { dimensions:

lon = 50 ; lat = 50 ;

variables:

byte Band1(lat, lon) ;

Band1:long_name = "GDAL Band Number 1" ; Band1:_Unsigned = "true" ; Band1:valid_range = 0s, 255s ; Band1:_FillValue = 0b ; Band1:grid_mapping = "crs" ;

char crs ;

crs:grid_mapping_name = "latitude_longitude" ; crs:longitude_of_prime_meridian = 0. ; crs:semi_major_axis = 6378137. ; crs:inverse_flattening = 298.257223563 ; crs:spatial_ref = "GEOGCS[\"unnamed ellipse\",DATUM[\"unknown\",SPHEROID[\"unnamed\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]]" ; crs:GeoTransform = "27 0.08 0 72 0 -0.04 " ;

double lat(lat) ;

lat:standard_name = "latitude" ; lat:long_name = "latitude" ; lat:units = "degrees_north" ;

double lon(lon) ;

lon:standard_name = "longitude" ; lon:long_name = "longitude" ; lon:units = "degrees_east" ;

global attributes:

:Conventions = "CF-1.5" ; :GDAL = "GDAL 2.0dev, released 2011/12/29" ; :history = "Tue Nov 06 10:07:19 2012: GDAL Create( gcps_pro_yes.nc, ... )" ;

data:

lat = 70.02, 70.06, 70.1, 70.14, 70.18, 70.22, 70.26, 70.3, 70.34, 70.38,

70.42, 70.46, 70.5, 70.54, 70.58, 70.62, 70.66, 70.7, 70.74, 70.78, 70.82, 70.86, 70.9, 70.94, 70.98, 71.02, 71.06, 71.1, 71.14, 71.18, 71.22, 71.26, 71.3, 71.34, 71.38, 71.42, 71.46, 71.5, 71.54, 71.58, 71.62, 71.66, 71.7, 71.74, 71.78, 71.82, 71.86, 71.9, 71.94, 71.98 ;

lon = 27.04, 27.12, 27.2, 27.28, 27.36, 27.44, 27.52, 27.6, 27.68, 27.76,

27.84, 27.92, 28, 28.08, 28.16, 28.24, 28.32, 28.4, 28.48, 28.56, 28.64, 28.72, 28.8, 28.88, 28.96, 29.04, 29.12, 29.2, 29.28, 29.36, 29.44, 29.52, 29.6, 29.68, 29.76, 29.84, 29.92, 30, 30.08, 30.16, 30.24, 30.32, 30.4, 30.48, 30.56, 30.64, 30.72, 30.8, 30.88, 30.96 ;

}

[etourigny]

I fail to see what is the problem.

The point in using WRITE_BOTTOM_UP=YES is to revert the y-axis in the netcdf output file, which the driver does correctly, as shown by ncdump output and also cdo.

$ cdo sinfo gcps_pro_yes.nc 

   File format: netCDF
    -1 : Institut Source   Param       Time Typ  Grid Size Num  Levels Num
     1 : unknown  unknown  -1          con  U8       2500   1       1   1
   Horizontal grids :
     1 : lonlat       > size      : dim = 2500  nlon = 50  nlat = 50
                        lon       : first = 27.04  last = 30.96  inc = 0.08  degrees_east
                        lat       : first = 70.02  last = 71.98  inc = 0.04  degrees_north
   Vertical grids :
     1 : surface                  : 0 
cdo sinfo: Processed 1 variable. ( 0.01s )

$ cdo sinfo gcps_pro_no.nc 

   File format: netCDF
    -1 : Institut Source   Param       Time Typ  Grid Size Num  Levels Num
     1 : unknown  unknown  -1          con  U8       2500   1       1   1
   Horizontal grids :
     1 : lonlat       > size      : dim = 2500  nlon = 50  nlat = 50
                        lon       : first = 27.04  last = 30.96  inc = 0.08  degrees_east
                        lat       : first = 71.98  last = 70.02  inc = -0.04  degrees_north
   Vertical grids :
     1 : surface                  : 0 
cdo sinfo: Processed 1 variable. ( 0.00s )

gdalinfo reports the same origin and (negative) y pixel size, because the GDAL model is north-up by default (origin is top-left corner, and y is in the south, or negative, direction). In fact, the netcdf driver (used by gdalinfo) reads the data in backwards when it detects a "bottom-up" file. So that is why you see the same origin (top-left) and negative pixel values.

The WRITE_BOTTOM_UP=YES creation option is not optimal for gdal use (because it will be slightly inefficient), but for other software which needs that to operate more efficiently (I can't name one right now). For example, ncview shows both files fine.

BTW why are you using that option?

Closing this bug as invalid.

[korosov]

The driver correctly writes data matrices and dimension variables.

It writes METADATA incorrectly. In a NetCDF file projection info is kept in attributes of one char variable usually called 'crs', and reference from each band by 'grid_mapping' attribute. This variable has attribute 'GeoTransform'. It is written incorrectly ("27 0.08 0 72 0 -0.04 ") in case of 'WRITE_BOTTOMUP=YES'. The last term (-0.04) should be positive because values of lat are increasing.

I noticed this problem when I tried to retrieve GeoTransform not using gdal (which calculates it from the dimension variable) but reading from metadata. That was needed for update of the DODS driver.

Now I will not rely on the 'GeoTransform' attribute in the 'crs' variable but will calculate it from dimension variable so priority for me is low.

[etourigny]

I see what you mean - GeoTransform and CF dimension variables are in contradiction.

I did not understand that you were referring to the netcdf metadata crs:GeoTransform , instead of the gdalinfo geotransform output.

BTW I didn't know that other software uses the crs GeoTransform metadata - where can I find info on this? I thought that that DODS would use standard netcdf coordinates variables.

Any application should look at the coordinate variables, not spatial_ref or GeoTransform attributes (which are not part of the CF/COARDS standard).

It is OK that gdalinfo reports negative y pixels, but not OK to report that in the CF metadata.

It is added merely as a convenience for gdal, not to be used outside of it. In fact it was added quite some time ago and should be removed, but I kept it for convenience and historical reasons. I should probably remove it in 2.0, but was waiting on some changes to the CF standard for better representing CRS information (like towgs84 params and such) inside CF attributes.

Another solution would be to re-calculate the GeoTransform origin and invert the y-axis when WRITE_BOTTOM_UP is used, but I don't have the time nor inclination to write the code and test.

I'm open to suggestions/patches on how to solve this.

[korosov]

We are updating the DODS driver and tried to read the 'GeoTransform' (if available in metadata) rather than calculate it. I don't know if other software uses it. I guess it can be just deleted from the attributes.

BTW, the DODS driver has a lot of common functionality with the netcdf. Ideally DODS should inherit from netcdf with ability to use libdap to access remote data and metadata. I hope we soon can come up with a patch and, maybe, rfc for the dods driver improvement. Opendap datasets are becoming more and more used nowadays and the way dods driver is implemented inhibits GDAL competitiveness, I believe.

[etourigny]

ah now I understand. You mean DODS driver in gdal - I thought you meant actual DODS server...

I see now that it does share common traits with the netcdf driver (at least for metadata) - although I have done significant changes to it from 1.9 onwards, hence the divergence. I confess I never even knew the gdal DODS drivver existed, nor looked at it.

I'll be willing to apply any patches you send, provided you give some test cases, ideally in the form of python autotests.

I'm curious how the GDAL metadata gets passed to the DODS driver from a netcdf file (generated by gdal)... what's a typical workflow, perhaps with an example?

You might want to stay away from any gdal-specific metadata, as it won't work for standard netcdf files, so there really is no point in using it. Files generated with gdal-1.9 netcdf driver are CF-1.6 compliant, that is output files are CF-1.6, but the driver doesn't support everything that CF-1.6 mandates (things like cell methods, dimension variable compression, etc.). You can probably use some code in the netcdf driver to read dimensions etc.

Like you propose DODS driver could subclass from netcdf driver, and the same could be done for gmt driver I suppose (it's been un-maintained for some time).

In the meantime - what would be the best solution for this problem? Removing GeoTransform seems the best choice, because there is no real added value and has some problems. spatial_ref will remain for some time, just because there is currently no standard way to store it in CF attributes.

[korosov]

I agree - removing GeoTransform is the best choice.

[etourigny]

this would probably be addressed with changes in cf 1.7, see ​https://cf-pcmdi.llnl.gov/trac/ticket/80 "Add missing CF parameters to translate Coordinate Reference System properties to/from OGC Well-Known Text format"

[Jukka Rahkonen]

The link in the last comment is dead. Did the comment mean that fixing that ticker number 80 would resolve this GDAL ticket without doing any work on the GDAL side?

[Markus Neteler]

Replying to jratike80:

The link in the last comment is dead. Did the comment mean that fixing that ticker number 80 would resolve this GDAL ticket without doing any work on the GDAL side?

FWIW: ​https://web.archive.org/web/20140418013544/http://cf-pcmdi.llnl.gov/trac/ticket/80

[etourigny]

Replying to jratike80:

The link in the last comment is dead. Did the comment mean that fixing that ticker number 80 would resolve this GDAL ticket without doing any work on the GDAL side?

The correct link is ​http://cf-trac.llnl.gov/trac/ticket/80´

And fixing this in gdal would require adding new metadata that are mentioned in the ticket, and removing useless ones like GeoTransform. But I never bothered implementing this in gdal as it has yet to be part of the official CF Conventions, even though it was accepted after long discussions.

[Even Rouault]

This ticket has been automatically closed because Trac is no longer used for GDAL bug tracking, since the project has migrated to GitHub. If you believe this ticket is still valid, you may file it to ​https://github.com/OSGeo/gdal/issues if it is not already reported there.