NetCDF driver assumes 'x' 'y' to be names of dimension variables

[korosov]

Problem

If a dimension variable has a name other than 'x' or 'y' the 'units' attribute is not taken into account while setting projection from CF in netCDFDataset::SetProjectionFromVar().

Example

The sample file is available here: ​http://thredds.met.no/thredds/catalog/cryoclim/met.no/osisaf-nh/catalog.html?dataset=cryoclim/met.no/osisaf-nh/osisaf-nh_aggregated_ice_concentration_nh_polstere-100_197810010000.nc

It looks nice in Godiva2 but dimension names are 'xc' and 'yc' and their units are 'km'. gdal cannot get the units and interprets it as 'm'. That results in wrong calculation of lat/lon.

Solution

Replace hardcoded "x" and "y" with corresponding poDS->papszDimName. Patch is attached.

[korosov]

Patch

[etourigny]

Thanks for the fix.

After applying a slightly modified patch, proj4 string shows units=km:

$ gdalsrsinfo NETCDF:"osisaf-nh_aggregated_ice_concentration_nh_polstere-100_197810010000.nc":ice_conc_avg

PROJ.4 : '+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=km +no_defs '

OGC WKT :
PROJCS["unnamed",
    GEOGCS["WGS 84",
[...]
    PARAMETER["false_northing",0],
    UNIT["Meter",1000,
        AUTHORITY["EPSG","9001"]]]

I have also found a bug in importFromProj4(), when using this dataset. It seems that proj4 import ignores units=km

$ gdalsrsinfo '+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=km +no_defs '

PROJ.4 : '+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs '

After both fixes, dataset shows fine in qgis (master).

Fixed both issues in trunk (r24772) and 1.9 (r24773).

[Even Rouault]

Etienne,

I believe that for km, it should rater be : oSRS.SetLinearUnits( "Kilometer", 1000.0 );

and the associated EPSG code should be 9036 (according to data/unit_of_measure.csv)

This will not change the translation to proj4 but should make the WKT to look more correct.

so I'd suggest the following patch:

Index: frmts/netcdf/netcdfdataset.cpp
===================================================================
--- frmts/netcdf/netcdfdataset.cpp	(révision 24774)
+++ frmts/netcdf/netcdfdataset.cpp	(copie de travail)
@@ -2631,8 +2631,8 @@
                     oSRS.SetAuthority( "PROJCS|UNIT", "EPSG", 9001 );
                 }
                 else if ( EQUAL(pszUnits,"km") ) {
-                    oSRS.SetLinearUnits( SRS_UL_METER, 1000.0 );
-                    oSRS.SetAuthority( "PROJCS|UNIT", "EPSG", 9001 );
+                    oSRS.SetLinearUnits( "Kilometer", 1000.0 );
+                    oSRS.SetAuthority( "PROJCS|UNIT", "EPSG", 9036 );
                 }
                 /* TODO check for other values */
                 // else 

[etourigny]

Hi Even,

shouldn't it also be kilometer in importFromProj4() ?

Also, there are conflicting definitions for m, which is correct? :

in ogr_srs_api.h :
#define SRS_UL_METER            "Meter"

in data/unit_of_measure.csv :
9001,metre,length,9001,1,1,Also known as International metre. SI standard unit.,ISO 1000.,OGP,1995/06/02,,0

[Even Rouault]

About "shouldn't it also be kilometer in importFromProj4() ? ", yes I missed that.

About meter vs metre, no idea where the reference should be found, if there's a reference actually ... So I wouldn't change anything for meter.

[etourigny]

fixed km UNITS name and authority in 1.9 (r24775) and trunk (24776).

EPSG 9036 is "kilometre" and EPSG 9001 is "metre".