NetCDF created by GDAL not CF-1.0 compliant for geographical grids

[bokhorst]

The NetCDF created by GDAL for a (non-projected) geographical grid does not seem to be CF-1.0 compliant, although the 'Conventions' attribute states otherwise.

The part that is not compliant is the grid_mapping_name "Geographics Coordinate System" that is not a valid grid_mapping_name in CF-1.0.

Here is an example of a GDAL-generated NetCDF (in CDL, only the header):

netcdf gdal_netcdf {
dimensions:
        x = 10 ;
        y = 10 ;
variables:
        char GDAL_Geographics ;
                GDAL_Geographics:Northernmost_Northing = 55. ;
                GDAL_Geographics:Southernmost_Northing = 60. ;
                GDAL_Geographics:Easternmost_Easting = 5. ;
                GDAL_Geographics:Westernmost_Easting = 0. ;
                GDAL_Geographics:spatial_ref = "GEOGCS[\"WGS
84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS
84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]"
;
                GDAL_Geographics:GeoTransform = "0 0.5 -0 55 0 0.5 " ;
                GDAL_Geographics:grid_mapping_name = "Geographics
Coordinate System" ;
                GDAL_Geographics:long_name = "Grid_latitude" ;
        float Band1(y, x) ;
                Band1:_FillValue = 0.f ;
                Band1:grid_mapping = "GDAL_Geographics" ;
                Band1:long_name = "GDAL Band Number 1" ;

// global attributes:
                :Conventions = "CF-1.0" ;
data:
}

[bokhorst]

Hi, I have created a patch that I'd like to share. Instead of using a grid mapping I now write out the lat/lon dimensions as data in the case of a geographical grid. Note that I center the data points in the grid cells instead of on its boundary.

[dnadeau]

Hi bokhorst,

There was an issue on writing lat/lon dimension when I created the driver. Some array can be extremely big and the conversion become time consuming.

I used the ESRI convention to specified the boundary to solve this problem for some users, but many software are now reading the lat/lon arrays you have created, and it might be worth doing it.

Thanks for your code, Denis

[etiennesky]

I see that this patch as not been included, could this be included in the 1.7 trunk?

What's stopping this feature from being included? As is, the driver is not CF-1.0 compliant.

My simple tests show that NetCDF->Geotiff correctly exports geographical information (included in the netCDF lat and lon information, but Geotiff->netCDF does not.

I'm willing to help dev and testing as I need accurate GeoTiff<->GDAL<->netcdf conversions.

Thanks Etienne

[etiennesky]

OK I've created a patch for 1.7.2 which is similar to the previous, with the following differences:

1 - GDAL_Geographics attributes are written, in addition to the lon and lat attributes 2 - "x" and "y" dimensions are named "lon" and "lat" 3 - "No Data" is set to NC_FILL_BYTE , because if not, then all 0 values are set to missing 4 - added support for all non-projected grids 2

---

I have tested this with one netCDF file and one GTIFF file, and a series of conversions shows that the data is intact and can be viewed in GrADS (a picky program when it comes to COARDS standards) and qview.

---

Change 4 : "if( ! oSRS.IsProjected() )" instead of "if( oSRS.IsGeographic() )" A quick fix to avoid code duplication, when there is no "GEOGCS" attribute This appeared necessary when doing netCDF->GTIFF->netCDF conversion for testing, as netCDF->GTIFF conversion does not add GEOGCS attributes.

As I don't have a file with a Projected grid, I have not implemented it for that type.

There is a lot of code duplication for the projected and geographic cases, which could be done outside of the 2 if statements.

[Kyle Shannon]

Although I haven't looked at this issue much, I do see that Denis was worried about substantial performance losses when creating and writing large lat/lon arrays. It does look like he thought this was a good idea for Geographic CS. I downloaded your patch and would like to test it on different datasets. Can you attach a small dataset for testing? thanks. kyle

[Kyle Shannon]

That is I applied your patch to my working copy and would like to test it...

[etiennesky]

Hmm not sure about a "small" dataset - the ones I have been using are GTIFF 24MB,5900X4200, 1km resolution.

No problem importing the file, a few seconds on my laptop. The lat/lon calculations and storage adds about a half-second time processing with my test file (LC15_amazon_biomass.tif).

You can probably have access through the official ftp server, if not I can compress it to about 3MB and attach it here.

​ftp://lbaworking.daac.ornl.gov/lba/land_use_land_cover_change/LC15_AGLB_Distribution_Map/data/LC15_amazon_biomass.tif ​ftp://lbaworking.daac.ornl.gov/lba/land_use_land_cover_change/LC15_AGLB_Distribution_Map/data/LC15_amazon_biomass.tif

[etiennesky]

I've made a few more modifications, and will upload the updated patch file (replaces previous one).

In addition to the first 3 changes in the previous patch:

1) only works on "Geographic" grids i.e. put back the if( oSRS.IsGeographic() ) statement
2) does not create the "GDAL_Geographics" variable, because this brakes the CF-1.0 standard, instead all values are set to global attributes "GDAL_Geographics::..."
3) invert latitude to conform with the "bottom-up" convention of netCDF. This has NOT been tested on datasets other than BYTE type, and should be extended to all coordinate types.
4) use double precision instead of single for the lat/lon values. This is needed for an accurate conversion of the grid data.

Can you test it with the above file and provide with a "Projected" file I could test it? How about non-geographic and non-projected? Do those files contain useful geo information that can be used to compute lat/lon values?

I think that most of the code in the 2 if sections ( oSRS.IsGeographic() and oSRS.IsProjected() ) should be merged for clarity and ease of use.

Sorry I'm not experienced with GDAL/GIS, coming from the Atmospheric Science, but I'm learning everyday. Thanks!

[etiennesky]

Another dataset which works (ARC/INFO format), but the file is large, is available at the following url:

​http://cdiac.ornl.gov/ftp/global_carbon/datasets.zip files datasets/c_10m/w001001.adf and datasets/c_5m/w001001.adf

attached is a tar.bz2 of the c_10m directory

This file is of FLOAT type, so it works with at least BYTE and FLOAT variables.

[etiennesky]

I've made a patch for recent version 1.8.0

Please test it and merge it with dev version!! This is an important issue, and the fix is quite simple. I've attached a test file from MODIS data, please test with it.

Also, I would like to see for myself the supposed performance problem with large lat/lon datasets. It really should not be that bad, considering there are far less lat and lon values than actual data in the Band(s). If someone can provide an example file I'd be happy to test it.

[etiennesky]

patch for 1.8.0

[etiennesky]

test file with modis data, in GTIFF format

[etiennesky]

I've made some small modifications to enhance metadata handling in netcdf import/export.

Also a fix for setting the CRS to WGS84 when importing a generic netcdf file.

File patch2-netcdf.txt contains these fixes, and file patch3-netcdf.txt also contains the previous fixes.

[Kyle Shannon]

I will take a look at the patch on monday, but I have a few quick comments/questions. Why should we assume WGS84 n generic import? Is it documented anywhere that CF-1.x conventions assume that datum? Just to be certain, pathc3-netcdf.txt contains all of the patches?

[etiennesky]

Thanks for taking the time.

As far as I understand, CF-1.x does not define datums, so as I understand it the WGS84 datum is always implicitly assumed. If you look at present code for projected CRS, WGS84 is used by default. I fixed the driver for setting WGS84 for regular lat/lon grids created with other software (i.e. no grid_mapping metadata). Datums are not considered in these cases, so I believe it is safe to assume WGS84. I've attached a small netcdf file as an example.

Yes, the 3rd patch includes all my fixes.

[etiennesky]

fixed in patch4-netcdf.txt improper handling of multi-band dataset.

If there are multiple bands with a NETCDF_VARNAME metadata, set the var names to <NETCDF_VARNAME><band_number>

[JoshVote]

I've just applied etiennesky's patch 'patch4-netcdf.txt' to revision 22366 of the gdal trunk with great results.

I've been using it to translate and ER mapper dataset into a NetCDF dataset for usage with the THREDDS WCS plugin.

Are there any plans for this fix to make it to the trunk?

[etiennesky]

I have ported my changes to gdal-1.8.1 with some differences.

There are 2 sets of fixes that I have tried to separate into 2 patches.

1) Fix problems related to raster export:

• add information on lat/lon for geographic grids, and remove custom gdal tags
• addresses the issue reported in ticket #4200 (coordinates shifting)
• fix an "upside-down" problem in exporting netcdf files (both geographical and projected) to support reading with other software
• when reading a netcdf file, check if it is "bottom first" (by checking the lat/lon for geographical grids, and a gdal tag "GDAL_Bottom_Up" for projected grids), maintaining backwards-compatibility

2) Make changes to the metadata handling (for proper import/export)

• Fix duplicate metatada as in issue #4204
• add variable data from a netcdf to Band metadata
• fix the variable name when exporting to netcdf

[etiennesky]

patch containing all fixes

[etiennesky]

patch for metadata

[etiennesky]

One issue which remains pending is how to deal with empty datum for lat/lon grids from other software. CF conventions leaves room for that in the spec (at least the geometric values), but no one uses it. Therefore most netcdf files are neither projected nor geographic (for GDAL).

I have modified the code to treat non-projected data as geographic to circumvent this problem. Another solution is to set WGS84 by default (when no other information is available), as I have suggested previously, however I am not sure this is valid.

Regards, Etienne

[Kyle Shannon]

I think it would be a good idea to apply the patch in 2 commits, one for the precision issues and another for the metadata issue. I will try to update all tickets involved with these issues. Thanks for the patch.

[etiennesky]

Kyle,

I will modify the first patch shortly, because the part that inverts the y-axis for projected CRS could break back-wards compatibility, and needs to be studied further. Even wrote on the mailing list that such changes require an RFC so I don't feel confident with that patch.

[etiennesky]

patch for metatada fixes

[etiennesky]

Added independent patches against svn trunk. You should be able to apply them sequentially. Sorry for multiple postings, but this one will work more easily.

1) patch-geo2.txt : Fixes export of geographical CRS

• adds lat/lon variables
• creates CF-1.5 compliant metadata for geographical crs (longitude_of_prime_meridian, semi_major_axis,inverse_flattening)
• y axis is inverted to be compatible with other software which expects "bottom-up"
• does not change projected CRS output (this will have to be addressed later)

2) patch-meta2.txt : Make changes to metadata handling (for proper import/export)

• Fix duplicate metatada as in issue #4204
• add metadata from a netcdf variable inside the corresponding GDAL Band instead of Global
• keep variable name when exporting to netcdf (after previous import)

[etiennesky]

same patch but without changes to dfNodata for NC_BYTE , just print a warning

[etourigny]

Added autotests to check for this bug in trunk (r23091), see bug #4233.

[etourigny]

Fixed in trunk (r23093).

Implemented changes in patch-geo2.txt with minor code cleanup. Autotests 18-20 now pass without problems, showing the fix is valid.

Conventions metadata has been updated to CF-1.2 because grid_mapping_name = "latitude_longitude" was introduced in CF-1.2.

An issue related to this fix (wiki:NetCDF_Improvements#Datumissues) - what should we do with any extra projection and datum information not supported by CF Conventions (EPSG, WKT, etc.)? For now they are dropped from files with geographic CRS, which is perfectly valid for this issue.

JoshVote and pds : can you update your svn and confirm that this fix works for you?

[hsumanto]

etienne: I've updated the svn and tested this fix work when converting dataset in epsg:4326 from geotiff -> netcdf and back from netcdf -> geotiff. The netcdf file can be viewed correctly with netcdf software.

This is the gdalinfo of the netcdf file.

Driver: netCDF/Network Common Data Format
Files: ep4326.nc
Size is 8990, 6284
Coordinate System is:
GEOGCS["unknown",
    DATUM["unknown",
        SPHEROID["Spheroid",6378137,298.257223563]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]
Origin = (144.496170150548124,-36.486543928430414)
Pixel Size = (0.000310407241244,-0.000310407241244)
Metadata:
  Band1#grid_mapping=crs
  Band1#long_name=GDAL Band Number 1
  crs#grid_mapping_name=latitude_longitude
  crs#inverse_flattening=298.257223563
  crs#longitude_of_prime_meridian=0
  crs#semi_major_axis=6378137
  lat#long_name=latitude
  lat#standard_name=latitude
  lat#units=degrees_north
  lon#long_name=longitude
  lon#standard_name=longitude
  lon#units=degrees_east
  NC_GLOBAL#Conventions=CF-1.2
  NC_GLOBAL#GDAL=GDAL 1.9dev, released 2011/01/18
  NC_GLOBAL#GDAL_AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  ( 144.4961702, -36.4865439) (144d29'46.21"E, 36d29'11.56"S)
Lower Left  ( 144.4961702, -38.4371430) (144d29'46.21"E, 38d26'13.71"S)
Upper Right ( 147.2867312, -36.4865439) (147d17'12.23"E, 36d29'11.56"S)
Lower Right ( 147.2867312, -38.4371430) (147d17'12.23"E, 38d26'13.71"S)
Center      ( 145.8914507, -37.4618435) (145d53'29.22"E, 37d27'42.64"S)
Band 1 Block=8990x1 Type=Byte, ColorInterp=Undefined
  NoData Value=0
  Metadata:
    grid_mapping=crs
    long_name=GDAL Band Number 1
    NETCDF_VARNAME=Band1

[etourigny]

r23615: fix upside-down export and import of grids without projection and geotransform (#2129, #4284)

[etourigny]

r23616 : only write grid_mapping attributes when bWriteGridMapping == TRUE (#2129, #2893 and #4294)

The problem was that files without grid_mapping had the spheroid parameters written as global metadata, which is incorrect. Also skip writing any projection info if ( !bIsProjected && !bWriteLonLat ).

[etourigny]

r23617 : partial revert of r23615 and set WRITE_BOTTOMUP default to YES for export and import of files without projection and GT (#2129, #4284)

previous fix broke many autotests, solution is to set bottom-up to TRUE by default, hen projection and GT are detected then the proper order is set (import and export). Hopefully this resolves a thorny question!