Opened 10 years ago

Closed 6 years ago

Last modified 6 years ago

#2129 closed task (fixed)

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

Reported by: bokhorst Owned by: dnadeau
Priority: normal Milestone:
Component: GDAL_Raster Version: 1.8.1
Severity: normal Keywords: netcdf, cf
Cc: Kyle Shannon, etiennesky, JoshVote, pds

Description

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:
}

Attachments (14)

netcdfdataset-1.5.0.patch (5.6 KB) - added by bokhorst 10 years ago.
netcdfdataset-1.7.2.patch (6.3 KB) - added by etiennesky 7 years ago.
netcdfdataset-1.7.2.patch2 (9.7 KB) - added by etiennesky 7 years ago.
patch1-netcdf.txt (15.3 KB) - added by etiennesky 6 years ago.
patch for 1.8.0
mcd2.tif (537.2 KB) - added by etiennesky 6 years ago.
test file with modis data, in GTIFF format
patch2-netcdf.txt (3.4 KB) - added by etiennesky 6 years ago.
patch3-netcdf.txt (18.0 KB) - added by etiennesky 6 years ago.
trmm1.nc (234.3 KB) - added by etiennesky 6 years ago.
patch4-netcdf.txt (20.4 KB) - added by etiennesky 6 years ago.
mod-all.txt (29.7 KB) - added by etiennesky 6 years ago.
patch containing all fixes
mod-1-out.txt (19.9 KB) - added by etiennesky 6 years ago.
mod-2-meta.txt (9.1 KB) - added by etiennesky 6 years ago.
patch for metadata
patch-meta2.txt (9.1 KB) - added by etiennesky 6 years ago.
patch for metatada fixes
patch-geo2.txt (17.2 KB) - added by etiennesky 6 years ago.
same patch but without changes to dfNodata for NC_BYTE , just print a warning

Download all attachments as: .zip

Change History (44)

comment:1 Changed 10 years ago by bokhorst

Type: defecttask
Version: svn-trunk1.5.0

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.

Changed 10 years ago by bokhorst

Attachment: netcdfdataset-1.5.0.patch added

comment:2 Changed 9 years ago by dnadeau

Owner: changed from warmerdam to dnadeau
Status: newassigned

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

comment:3 Changed 7 years ago by Kyle Shannon

Keywords: netcdf cf cf-1.0 added; NetCDF CF CF-1.0 removed

comment:4 Changed 7 years ago by 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

comment:5 Changed 7 years ago by etiennesky

Version: 1.5.01.7.2

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.

Changed 7 years ago by etiennesky

Attachment: netcdfdataset-1.7.2.patch added

comment:6 Changed 7 years ago by 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

comment:7 Changed 7 years ago by Kyle Shannon

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

comment:8 Changed 7 years ago by 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

comment:9 Changed 7 years ago by 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!

Changed 7 years ago by etiennesky

Attachment: netcdfdataset-1.7.2.patch2 added

comment:10 Changed 7 years ago by 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.

comment:11 Changed 6 years ago by etiennesky

Version: 1.7.21.8.0

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.

Changed 6 years ago by etiennesky

Attachment: patch1-netcdf.txt added

patch for 1.8.0

Changed 6 years ago by etiennesky

Attachment: mcd2.tif added

test file with modis data, in GTIFF format

comment:12 Changed 6 years ago by 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.

Changed 6 years ago by etiennesky

Attachment: patch2-netcdf.txt added

Changed 6 years ago by etiennesky

Attachment: patch3-netcdf.txt added

comment:13 Changed 6 years ago by 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?

comment:14 Changed 6 years ago by 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.

Changed 6 years ago by etiennesky

Attachment: trmm1.nc added

Changed 6 years ago by etiennesky

Attachment: patch4-netcdf.txt added

comment:15 Changed 6 years ago by 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>

comment:16 Changed 6 years ago by 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?

comment:17 Changed 6 years ago by Even Rouault

Cc: Kyle Shannon added

comment:18 Changed 6 years ago by 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

Changed 6 years ago by etiennesky

Attachment: mod-all.txt added

patch containing all fixes

Changed 6 years ago by etiennesky

Attachment: mod-1-out.txt added

Changed 6 years ago by etiennesky

Attachment: mod-2-meta.txt added

patch for metadata

comment:19 Changed 6 years ago by etiennesky

Version: 1.8.01.8.1

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

comment:20 Changed 6 years ago by 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.

comment:21 Changed 6 years ago by 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.

Changed 6 years ago by etiennesky

Attachment: patch-meta2.txt added

patch for metatada fixes

comment:22 Changed 6 years ago by 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)

Changed 6 years ago by etiennesky

Attachment: patch-geo2.txt added

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

comment:23 Changed 6 years ago by etiennesky

Cc: etiennesky added

comment:24 Changed 6 years ago by etourigny

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

comment:25 Changed 6 years ago by etourigny

Cc: JoshVote pds added
Keywords: cf-1.0 removed

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?

comment:26 Changed 6 years ago by 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

comment:27 Changed 6 years ago by etourigny

Resolution: fixed
Status: assignedclosed

comment:28 Changed 6 years ago by etourigny

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

comment:29 Changed 6 years ago by 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 ).

comment:30 Changed 6 years ago by 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!

Note: See TracTickets for help on using tickets.