Opened 9 years ago

Closed 9 years ago

#4769 closed defect (fixed)

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

Reported by: korosov Owned by: warmerdam
Priority: normal Milestone: 1.9.2
Component: default Version: svn-trunk
Severity: normal Keywords: NetCDF, dimensions
Cc: etourigny

Description

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.

Attachments (1)

4769_netcdf_dimensions.diff (869 bytes) - added by korosov 9 years ago.
Patch

Download all attachments as: .zip

Change History (7)

Changed 9 years ago by korosov

Attachment: 4769_netcdf_dimensions.diff added

Patch

comment:1 Changed 9 years ago by Even Rouault

Cc: etourigny added

comment:2 Changed 9 years ago by etourigny

Resolution: fixed
Status: newclosed

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).

comment:3 Changed 9 years ago by Even Rouault

Resolution: fixed
Status: closedreopened

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 

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

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

comment:6 Changed 9 years ago by etourigny

Resolution: fixed
Status: reopenedclosed

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

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

Note: See TracTickets for help on using tickets.