netCDF: implement Get/SetUnitType() using the standard "units" attribute

autotest/gdrivers/netcdf.py

@@ -2128,6 +2128,29 @@ def netcdf_56():
 
     return 'success'
 
+###############################################################################
+#check for UnitType set/get.
+def netcdf_57():
+
+    if gdaltest.netcdf_drv is None:
+        return 'skip'
+
+    # get
+    ds = gdal.Open( 'data/unittype.nc' )
+
+    unit = ds.GetRasterBand( 1 ).GetUnitType();
+
+    if unit != 'm/s':
+        gdaltest.post_reason( 'Incorrect unit(%s)' % unit )
+        return 'fail'
+
+    ds = None
+
+    # set
+    tst = gdaltest.GDALTest( 'NetCDF', 'unittype.nc', 1, 4672 )
+
+    return tst.testSetUnitType()
+
 ###############################################################################
 
 ###############################################################################
@@ -2189,7 +2212,8 @@ def netcdf_56():
     netcdf_53,
     netcdf_54,
     netcdf_55,
-    netcdf_56
+    netcdf_56,
+    netcdf_57
 ]
 
 ###############################################################################

autotest/pymod/gdaltest.py

@@ -985,6 +985,50 @@ def testSetDescription(self):
 
         return 'success'
 
+    def testSetUnitType(self):
+        if self.testDriver() == 'fail':
+            return 'skip'
+
+        src_ds = gdal.Open( 'data/' + self.filename )
+        xsize = src_ds.RasterXSize
+        ysize = src_ds.RasterYSize
+
+        new_filename = 'tmp/' + self.filename + '.tst'
+        new_ds = self.driver.Create( new_filename, xsize, ysize, 1,
+                                     src_ds.GetRasterBand(self.band).DataType,
+                                     options = self.options  )
+        if new_ds is None:
+            post_reason( 'Failed to create test file using Create method.' )
+            return 'fail'
+
+        unit = 'mg/m3'
+        if new_ds.GetRasterBand(1).SetUnitType( unit ) is not gdal.CE_None:
+            post_reason( 'Failed to set unit type.' )
+            return 'fail'
+
+        src_ds = None
+        new_ds = None
+
+        new_ds = gdal.Open( new_filename )
+        if new_ds is None:
+            post_reason( 'Failed to open dataset: ' + new_filename )
+            return 'fail'
+
+        new_unit = new_ds.GetRasterBand(1).GetUnitType()
+        if new_unit != unit:
+            print('')
+            print('old = ', unit)
+            print('new = ', new_unit)
+            post_reason( 'Did not get expected unit type.' )
+            return 'fail'
+
+        new_ds = None
+
+        if gdal.GetConfigOption( 'CPL_DEBUG', 'OFF' ) != 'ON':
+            self.driver.Delete( new_filename )
+
+        return 'success'
+
 
 def approx_equal( a, b ):
     a = float(a)

gdal/frmts/netcdf/netcdfdataset.cpp

@@ -113,6 +113,7 @@ class netCDFRasterBand : public GDALPamRasterBand
     double      adfValidRange[2];
     double      dfScale;
     double      dfOffset;
+    CPLString   osUnitType;
     bool        bSignedData;
     bool        bCheckLongitude;
 
@@ -154,6 +155,8 @@ class netCDFRasterBand : public GDALPamRasterBand
     virtual CPLErr SetOffset( double );
     virtual double GetScale( int * );
     virtual CPLErr SetScale( double );
+    virtual const char *GetUnitType();
+    virtual CPLErr SetUnitType( const char * );
     virtual CPLErr IReadBlock( int, int, void * );
     virtual CPLErr IWriteBlock( int, int, void * );
 };
@@ -430,6 +433,11 @@ netCDFRasterBand::netCDFRasterBand( netCDFDataset *poNCDFDS,
         CPLTestBool(CPLGetConfigOption("GDAL_NETCDF_CENTERLONG_180", "YES"))
         && NCDFIsVarLongitude( cdfid, nZId, NULL );
 
+/* -------------------------------------------------------------------- */
+/* Attempt to fetch the units attribute for the variable and set it.    */
+/* -------------------------------------------------------------------- */
+    SetUnitType( GetMetadataItem( CF_UNITS ) );
+
 /* -------------------------------------------------------------------- */
 /*      Check for variable chunking (netcdf-4 only)                     */
 /*      GDAL block size should be set to hdf5 chunk size                */
@@ -747,6 +755,52 @@ CPLErr netCDFRasterBand::SetScale( double dfNewScale )
     return CE_None;
 }
 
+/************************************************************************/
+/*                            GetUnitType()                             */
+/************************************************************************/
+
+const char *netCDFRasterBand::GetUnitType()
+
+{
+    if( osUnitType.size() > 0 )
+        return osUnitType;
+    else
+        return GDALRasterBand::GetUnitType();
+}
+
+/************************************************************************/
+/*                           SetUnitType()                              */
+/************************************************************************/
+
+CPLErr netCDFRasterBand::SetUnitType( const char* pszNewValue )
+
+{
+    CPLMutexHolderD(&hNCMutex);
+
+    osUnitType = (pszNewValue != NULL ? pszNewValue : "");
+
+    if ( osUnitType.size() > 0 ) {
+        /* write value if in update mode */
+        if ( poDS->GetAccess() == GA_Update ) {
+
+            /* make sure we are in define mode */
+            ( ( netCDFDataset * ) poDS )->SetDefineMode( TRUE );
+
+            int status = nc_put_att_text( cdfid, nZId, CF_UNITS,
+                                          osUnitType.size(),
+                                          osUnitType.c_str() );
+
+            NCDF_ERR(status);
+            if ( status == NC_NOERR )
+                return CE_None;
+
+            return CE_Failure;
+        }
+    }
+
+    return CE_None;
+}
+
 /************************************************************************/
 /*                           GetNoDataValue()                           */
 /************************************************************************/