netcdf driver does not read and write CF-1.0 coordinate system properly for projected CRS

[warmerdam]

If I generate a netcdf file from a UTM GeoTIFF file like this:

gdal_translate -srcwin 0 0 20 20 utm.tif out.nc -of netcdf

and then remove the spatial_ref and geotransform metadata items which are non-standard GDAL items (using ncdump/vi/ncgen) I end up with what should have been a valid CF-1.0 file that GDAL does not recognise the coordinate system from, nor extract the geotransform.

Driver: netCDF/Network Common Data Format
Files: out.nc
       out.nc.aux.xml
Size is 20, 20
Coordinate System is `'
Metadata:
  NC_GLOBAL#Conventions=CF-1.0
  NC_GLOBAL#AREA_OR_POINT=Area
  Band1#grid_mapping=transverse_mercator
  Band1#long_name=GDAL Band Number 1
  transverse_mercator#Northernmost_Northing=3.75132e+06
  transverse_mercator#Southernmost_Northing=3.75012e+06
  transverse_mercator#Easternmost_Easting=441920
  transverse_mercator#Westernmost_Easting=440720
  transverse_mercator#grid_mapping_name=transverse_mercator
  transverse_mercator#latitude_of_projection_origin=0.000000e+00
  transverse_mercator#central_meridian=-1.170000e+02
  transverse_mercator#scale_factor_at_central_meridian=9.996000e-01
  transverse_mercator#false_easting=5.000000e+05
  transverse_mercator#false_northing=0.000000e+00
  AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,   20.0)
Upper Right (   20.0,    0.0)
Lower Right (   20.0,   20.0)
Center      (   10.0,   10.0)
Band 1 Block=20x1 Type=Byte, ColorInterp=Gray
  NoData Value=0
  Metadata:
    NETCDF_VARNAME=Band1

We need to review the GDAL NetCDF driver logic for reading this CF-1.0 info and try and correct it. It is also possible we are writing something wrong.

[kyle]

gdal does write something wrong, kind of. GDAL does not write CF-1.0 tags. Instead it writes its own tags, the easting/northing tags, geotransform, and spatial_ref. These cover the values CF has for the grid_mapping. The driver should probably be updated to write CF as well as gdal. Another issue is that the driver may assume that if one gdal tag is written, it contains the others, and vice versa. Right now, mixing tags is not supported as far as I can see. In the example above the geotransform can be calculated from the easting/northing tags, but it isn't, the driver searches for geo_transform. This issue probably falls under the work I am trying to do to fix some tag issues with the driver. I hope to have a solid fix sometime in the near future.

[hsumanto]

tile.nc can not be recognized as grid by toolsUI Gridviewer

[hsumanto]

tile-fixed-new.nc can be recognized as grid by toolsUI Gridviewer after fixing projection coordinates x and y and projection standard parallel values

[hsumanto]

dataset info for tile.nc where it can not find x and y axis and also projected CRS

[hsumanto]

tile-fixed-new.nc dataset info

[hsumanto]

There has been a patch in ticket #2129 to fix geographic CRS: -adds lat/lon variables -creates CF-1.5 compliant metadata for geographical crs (longitude_of_prime_meridian, semi_major_axis,inverse_flattening)

but not for projected CRS.

I am currently having issue trying to view the NetCDF dataset (Note: this file was created by gdal_translate geotiff LANDSAT dataset into NetCDF and has projected CRS) using netCDF toolsUI Gridviewer as it is not recognized as a grid. See  http://www.unidata.ucar.edu/mailing_lists/archives/netcdf-java/2011/msg00123.html

Looking at this documentation ( http://www.unidata.ucar.edu/software/netcdf-java/reference/StandardCoordinateTransforms.html), in order for grid data to be recognized as a grid by Gridviewer, the NetCDF dataset must: a) define x and y projection coordinate variables, using the correct projection units on the projection plane. b) define your projection dummy variable which has an attribute "grid_mapping_name" c) refer to the projection in your data variables with the "grid_mapping" attribute.

Checking on the GDAL netCDF, it only define x and y dimensions but not x and y projection coordinate variables, thus, this is the first issue (point a is not followed).

I wondering whether other people have experienced the same issue as what I am having.

Then I just tried to: -generate its CDL by ncdump the GDAL netCDF file into CDL -and the define x and y projection coordinate variables (as shown below) into the CDL

double y(y) ;

y:units = "m" ; y:standard_name = "projection_y_coordinate" ; y:long_name = "y coordinate of projection" ;

double x(x) ;

x:units = "m" ; x:standard_name = "projection_x_coordinate" ; x:long_name = "x coordinate of projection" ;

-regenerate the netCDF file from the modified CDL using ncgen -populate the y coordinate variable in the NetCDF file with

ymin, ymin+resolution, ymin+2*resolution, ...

-populate the x coordinate variable with

xmin, xmin+resolution, xmin+2*resolution, ...

-if I then viewed the modified NetCDF dataset again using Gridviewer in NetCDF Tools (4.2), now it can be viewed correctly and is recognized as a grid.

Note: xmin, ymin are the lower left corner and resolution is the raster pixel size.

Thus, my questions: i)Is there a specific reason why the current driver do not define the x and y projection coordinate variables but only defining x and y dimensions. For the dataset to be recognized as a grid by tool such as netCDF toolsUI, we have to add projection coordinate variables and their correct data values into the netCDF file.

ii)For lambert conformal conic projection with more than one standard parallel values, GDAL netCDF file is showing lambert_conformal_conic:standard_parallel_1 = -18.f ; lambert_conformal_conic:standard_parallel_2 = -36.f ;

I suspected standard_parallel_1 and standard_parallel_2 are not recognized by NetCDF ToolsUI (but strangely the file does in deed pass CF compliant checker) because by changing it to lambert_conformal_conic:standard_parallel = -18.f,-36.f ;

I managed to get the netcdf file to be recognized as a grid by Gridviewer in NetCDF Tools (4.2).

I would like to attached both files below but unable to do so because of the size: tile.nc - GDAL netCDF file tile-fixed-new.nc - the one I just tried to fix and can be viewed in netCDF toolsUI

I did attach some screenshots and the dataset info I get using Gridviewer in NetCDF Tools (4.2).

[etiennesky]

Hendy,

can you add compressed (zip or gzip) versions of the gtiff and perhaps the netcdf files? Perhaps you can make a spatial subset to make them smaller, or upload them somewhere else.

Regarding the information you sent about the files, it is more standard to use "ncdump -h file.nc" as everybody has access to ncdump. Your output is somewhat different.

The CF conventions checker seemingly only checks for valid names and doesn't seem to look very deeply into grid structure, so I don't think it is a robust validator for the moment.

For your questions

i) the driver writes what it needs in order to be able to read the geotransform and projection info using custom tags, it needs to be updated.

ii) This has been discussed in ticket #3324 . It seems that import supports the 2 ways of defining standard_parallel, but output is limited to standard_parallel_1 and standard_parallel_2 (the GDAL way). The fix would be to output a single standard_parallel with 2 values (as the CF convention stipulates), which is already supported for input.

[etiennesky]

patch to create a single standard_parallel variable that is CF compliant

[etiennesky]

Hendy, can you please apply the patch-std-parallel.txt patch and see if it resolves your second issue?

[hsumanto]

Hi etiennesky, I have tried the patch for std parallel, it works OK. Thanks for that.

This is the ncdump result of the patch netcdf file.

netcdf melb3112-patch {
dimensions:
        x = 198 ;
        y = 99 ;
variables:
        char lambert_conformal_conic ;
                lambert_conformal_conic:Northernmost_Northing = -4218968.14605944 ;
                lambert_conformal_conic:Southernmost_Northing = -4221948.66130479 ;
                lambert_conformal_conic:Easternmost_Easting = 1060889.92068945 ;
                lambert_conformal_conic:Westernmost_Easting = 1054928.89019874 ;
                lambert_conformal_conic:spatial_ref = "PROJCS[\"GDA94 / Geoscience Australia Lambert\",GEOGCS[\"GDA94\",DATUM[\"Geocentric_Datum_of_Australia_1994\",SPHEROID[\"GRS 1980\",6378137,298.2572221010002,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"6283\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4283\"]],PROJECTION[\"Lambert_Conformal_Conic_2SP\"],PARAMETER[\"standard_parallel_1\",-18],PARAMETER[\"standard_parallel_2\",-36],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",134],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AUTHORITY[\"EPSG\",\"3112\"]]" ;
                lambert_conformal_conic:GeoTransform = "1.05493e+06 30.1062 0 -4.21897e+06 0 -30.1062 " ;
                lambert_conformal_conic:grid_mapping_name = "lambert_conformal_conic" ;
                lambert_conformal_conic:latitude_of_projection_origin = 0. ;
                lambert_conformal_conic:longitude_of_central_meridian = 134. ;
                lambert_conformal_conic:false_easting = 0. ;
                lambert_conformal_conic:false_northing = 0. ;
                lambert_conformal_conic:standard_parallel = -18., -36. ;
        byte Band1(y, x) ;
                Band1:grid_mapping = "lambert_conformal_conic" ;
                Band1:long_name = "GDAL Band Number 1" ;

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

I have also attached the geotiff file which can be used for geotiff -> netcdf conversion testing.

The only thing left is to add projection coordinate variables (x and y) and their values into the netcdf file.

I can could help you test the patch once it is done and I will also use this dataset to test whether it fixes the inverted y axis for projected CRS.

[hsumanto]

small geotif dataset which could be used for testing geotif -> netCDF

[hsumanto]

The netcdf file which I have fixed by adding projection coordinate variables x and y

[hsumanto]

I have added the netcdf which I have fixed myself by doing the steps I have described above.

I am also wondering why I am still experiencing the inverted y-axis issue mentioned in ticket #2129 even though I have applied patch-geo2.txt when viewing the netcdf file using NetCDF Tools (4.2). The dataset is using projected CRS.

[etiennesky]

Thanks for the files.

I did not fix the y-axis for projected CRS, because this topic is related to geographic crs. As the driver expects top-down for projected crs, I did not want to create an internal incompatibility. This is one of the issues to address in the fixes for projected crs.

If you want to flip the y-axis, just apply the following around line 369 (in the part that deals with projected crs. However, the resulting file will be upside-down in gdal and related software.

/* netcdf standard is bottom-up, but leave it top first for now */

- bBottomUp = FALSE; + bBottomUp = TRUE;

[etiennesky]

hmm that got mangled...

       /* netcdf standard is bottom-up, but leave it top first for now */
-       bBottomUp = FALSE;
+       bBottomUp = TRUE;

[hsumanto]

Thanks for that, with that fix the dataset now looks correct when viewed with netCDF ToolsUI.

Now just need to wait projection coordinate variables together with theirs values being added into netCDF file.

Would you be slowly adding all these fixes into trunk once you get the commit access? As there are quite a few good fixes, you had made so far.

[hsumanto]

I added the below into the wiki as it would be good to get it fixed for projected crs once everyone agrees on the total solution which still maintains the compatibility between GDAL tools and netCDF tools.

"fix the y-axis for projected CRS, because currently the driver (GDAL) expects top-down for projected crs but netCDF is bottom-up approach. This issue need to be addressed for projected crs so that the resulting file in geotiff -> netcdf conversion (and vice versa netcdf -> geotiff) will not be upside-down when viewed in both gdal and netcdf software" 

[pds]

Also re this issue: recall that we still need to address the issue of saving Datum information with projected CRS's in NetCDF, or else the referenced files will be displaced.

I know datum issues are listed at wiki:NetCDF_Improvements#Datumissues, but adding this as a reminder as it's a potential big 'gotcha' for doing things like Thredds WMS, even after applying the earlier patches. E.g. to reference with WGS84 as the datum you'd add:

    crs:longitude_of_prime_meridian = 0.0 ;
    crs:semi_major_axis = 6378137.0 ;
    crs:inverse_flattening = 298.257223563 ;

[etiennesky]

hsumanto: if I get commit access I will start with incorporating the compliance checker (ticket #4233) and do some autotests, to make sure that any fixes really make a difference in CF-compliance and don't break other things. Then I will incorporate the geographic CRS fixes and metadata fixes, as well as a few simple fixes (e.g. #3324 and #2379).

The projected CRS fixes will have to wait a "final" solution and code. Any contribution from you guys is welcome or course, and thanks for your help!

pds: added a reminder in the CF topic in the wiki. Is the fix for geographic CRS good enough for Thredds? If not what is missing?

Just wondering, could you help write the code for proper projected import and export once we are settled? I don't have sufficient experience with projections to write it on my own. Thanks

[hsumanto]

Cool, happy to know you will add them into trunk.

The projected CRS fixes will have to wait a "final" solution and code. Any >>contribution from you guys is welcome or course, and thanks for your help!

Most of issues which I have encountered so far have been added to projected crs issue list in wiki:NetCDF_Improvements. I will try to learn the current driver first.

Is the fix for geographic CRS good enough for Thredds? If not what is missing?

I have tried the fix for geographic CRS and when viewing the netcdf dataset in Thredds, it can be viewed without any problem.

[hsumanto]

etienne: Just wondering whether you will commit the patch-std-parallel.txt patch into trunk.

[etourigny]

I have to cleanup the std-parallel fix before submitting it and make sure it will be the final solution.

I have done some preliminary work on fully compliant export for projected crs, sent you the info by email.

[pds]

Another CF-1 Map projection found similar to the std_parallel issue for LCC: Albers equal area has a problem with the 'latitude_of_projection_origin' attrib.

Relevant GDALInfo section of the original sample Tiff file:

    PROJECTION["Albers_Conic_Equal_Area"],
    PARAMETER["standard_parallel_1",29.5],
    PARAMETER["standard_parallel_2",45.5],
    PARAMETER["latitude_of_center",23],
    PARAMETER["longitude_of_center",-96],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],

Once translated to NetCDF (including with the std_parallel patch applied), this section is the same based on the WKT, but the "latitude_of_center" failed to make it across to 'latitude_of_projection_origin' in CF-1, ncdump shown:

		albers_conical_equal_area:longitude_of_central_meridian = -96. ;
		albers_conical_equal_area:false_easting = 0. ;
		albers_conical_equal_area:false_northing = 0. ;
		albers_conical_equal_area:standard_parallel = 29.5, 45.5 ;

So I had to manually add the below to get it to read properly in ToolsUI/IDV:

        albers_conical_equal_area:latitude_of_projection_origin = 23. ;

[etourigny]

Pat,

assuming that GDAL latitude_of_center maps to CF latitude_of_projection_origin (which seems to be the case), you have to add a relevant definition in netcdfdataset.h

{"latitude_of_projection_origin", SRS_PP_LATITUDE_OF_ORIGIN },

+ {"latitude_of_projection_origin", SRS_PP_LATITUDE_OF_CENTER },

{FALSE_EASTING, SRS_PP_FALSE_EASTING },

The string values for the SRS_PP_* and SRS_PT* are given in ogr/ogr_srs_api.h

Please try that, thanks!

[etourigny]

actually it would be cleaner to use the following

    {SCALE_FACTOR, SRS_PP_SCALE_FACTOR },   
    {STD_PARALLEL_1, SRS_PP_STANDARD_PARALLEL_1 },
    {STD_PARALLEL_2, SRS_PP_STANDARD_PARALLEL_2 },
    {LONG_CENTRAL_MERIDIAN, SRS_PP_CENTRAL_MERIDIAN },
    {LONG_CENTRAL_MERIDIAN, SRS_PP_LONGITUDE_OF_CENTER },
    {LON_PROJ_ORIGIN, SRS_PP_LONGITUDE_OF_ORIGIN }, 
    {LAT_PROJ_ORIGIN, SRS_PP_LATITUDE_OF_ORIGIN }, 
    {LAT_PROJ_ORIGIN, SRS_PP_LATITUDE_OF_CENTER }, 

Please test this, it might not be completely accurate...

[pds]

Update: I've done some testing with Etienne's suggested changes above, and it shows:

Depending on how SRS_PP_LONGITUDE_OF_CENTER mapping is defined, we can now export with CF-1.5 compliance to following projections:

• Albers Equal Area:
• Azimuthal equidistant:
• Lambert azimuthal equal area
• Lambert conformal
• Lambert Cylindrical Equal Area
• Transverse Mercator

However, there are still parameter problems for several other projections, mainly to with how various longitude parameters are mapped:

• Mercator, Orthographic, Polar stereographic, Stereographic

I will attach an RST file showing details of this testing (I've been using Python scripts to test all projections).

This shows that the mapping of OGC parameter names GDAL uses, to CF-1 parameter names, is dependent per-projection used, so to me that means the big universal mapping array in source:trunk/gdal/frmts/netcdf/netcdfdataset.h (poNetcdfSRS) is not going to be sufficient.

Seems to me a workable approach would be a 2D array, listing the mapping of parameter name mappings per-projection. This could possibly help reduce duplicate logic in the SetProjection?() function in netcdfdataset.cpp, which currently requires a big if statement with different logic for each projection.

[pds]

Notes on mapping between projections, and current status of testing.

[etourigny]

changes for projections in trunk r23449: Update comments about Polar Stereographic projection r23450: export SCALE_FACTOR_ORIGIN for LCC-1SP so it can be read back into GDAL

more changes coming today

[etourigny]

see #3324 for changes in LCC-1SP

[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]

r23627 : support import of polar stereographic variant without standard parallel (#2893)

standard_parallel can be computed from scale_factor using Snyder eq. 22-7 , see the header file and wik:NetCDF_ProjectionTestingStatus for more details

[etourigny]

leaving this bg open for now, pending

1) resolution of issues with 3 projections (see  http://mailman.cgd.ucar.edu/pipermail/cf-metadata/2011/027605.html )

2) changes to CF standard regarding missing WKT parameters such as DATUM and TOWGS84 (see  http://mailman.cgd.ucar.edu/pipermail/cf-metadata/2011/027499.html ).

[etourigny]

2) 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"