netCDF: allow Reading from CF-1.8 Encoded Geometries (#1287) (#1541)

gdal/doc/source/drivers/vector/netcdf.rst

@@ -16,9 +16,12 @@ representing scientific data.
 
 The driver handles the "point" and "profile" `feature
 types <http://cfconventions.org/cf-conventions/v1.6.0/cf-conventions.html#_features_and_feature_types>`__
-of the CF 1.6 convention. It also supports a more custom approach for
+of the CF 1.6 convention. For CF-1.7 and below (as well as non-CF files), it also supports a more custom approach for
 non-point geometries.
 
+The driver also supports reading from CF-1.8 convention compliant files that
+have simple geometry information encoded within them.
+
 Driver capabilities
 -------------------
 
@@ -93,14 +96,19 @@ double (with units="seconds since 1970-1-1 0:0:0") DateTime
 
 Layers
 ~~~~~~
-
-Generally a single netCDF file is viewed as a single OGR layer, provided
-that it contains only mono-dimensional variables, indexed by the same
-dimension (or bi-dimensional variables of type char). For netCDF v4
-files with multiple groups, each group may be seen as a separate OGR
-layer.
-
-On writing, the `MULTIPLE_LAYERS <#MULTIPLE_LAYERS>`__ dataset creation
+In the CF-1.8 compliant driver, a single layer corresponds to a single
+**geometry container** within a CF-1.8 compliant netCDF file. A geometry container, per
+the CF-1.8 specification, is referred to by another variable
+(presumably a data variable) through the **geometry** attribute. When reading
+a CF-1.8 compliant netCDF file, all geometry containers within the netCDF file
+will be present in the opened dataset as separate layers.
+
+When working with files made with older versions of the driver (pre CF-1.8),
+a single netCDF file generally corresponds to a single OGR layer,
+provided that it contains only mono-dimensional variables,
+indexed by the same dimension (or bi-dimensional variables of type char).
+For netCDF v4 files with multiple groups, each group may be seen as a separate OGR
+layer. On writing, the `MULTIPLE_LAYERS <#MULTIPLE_LAYERS>`__ dataset creation
 option can be used to control whether multiple layers is disabled, or if
 multiple layers should go in separate files, or separate groups.
 
@@ -123,6 +131,13 @@ written as netCDF v4 variable length strings.
 
 Geometry
 ~~~~~~~~
+Supported feature types when reading from a CF-1.8 convention compliant netCDF file
+include OGRPoint, OGRLineString, OGRPolygon, OGRMultiPoint, OGRMultiLineString, and
+OGRMultiPolygon. Due to slight ambiguities present in the CF-1.8 convention concerning
+Polygons versus MultiPolygons, the driver will in most cases default to assuming a MultiPolygon
+for the geometry of a layer with **geometry_type** polygon. The one exception where a Polygon type
+will be used is when the attribute **part_node_count** is not present within that layer's geometry container.
+Per convention requirements, the driver supports reading from geometries with X, Y, and Z axes.
 
 Layers with a geometry type of Point or Point25D will cause the implicit
 creation of x,y(,z) variables for projected coordinate system, or
@@ -316,12 +331,23 @@ The following example shows all possibilities and precedence rules:
 The effect on the output can be checked by running the **ncdump**
 utility
 
-See Also:
----------
+Limitation(s)
+-------------
+Only one grid mapping may be associated with an entire dataset of netCDFLayers.
+Attempting use more than one grid mapping, will result in the last set grid mapping
+overwriting all others.
+
+Furthermore, CF-1.8 grid mappings are currently not supported. Support for CF-1.8 grid
+mappings will be introduced in the near future.
+
+Further Reading
+---------------
 
 -  :ref:`Raster side of the netCDF driver. <raster.netcdf>`
 -  `NetCDF CF-1.6
    convention <http://cfconventions.org/cf-conventions/v1.6.0/cf-conventions.html>`__
+-  `NetCDF CF-1.8
+   convention draft <https://github.com/cf-convention/cf-conventions/blob/master/ch07.adoc>`__
 -  `NetCDF compiled
    libraries <http://www.unidata.ucar.edu/downloads/netcdf/index.jsp>`__
 -  `NetCDF

gdal/frmts/netcdf/GNUmakefile

@@ -1,7 +1,7 @@
 
 include ../../GDALmake.opt
 
-OBJ	=	netcdfdataset.o netcdflayer.o netcdfwriterconfig.o gmtdataset.o
+OBJ	=	netcdfdataset.o netcdflayer.o netcdfwriterconfig.o gmtdataset.o netcdfsg.o netcdflayersg.o
 
 XTRA_OPT =
 ifeq ($(NETCDF_HAS_NC4),yes)

gdal/frmts/netcdf/makefile.vc

@@ -3,7 +3,7 @@ GDAL_ROOT	=	..\..
 
 !INCLUDE $(GDAL_ROOT)\nmake.opt
 
-OBJ		=	netcdfdataset.obj netcdflayer.obj netcdfwriterconfig.obj gmtdataset.obj
+OBJ		=	netcdfdataset.obj netcdflayer.obj netcdfwriterconfig.obj gmtdataset.obj netcdfsg.obj netcdflayersg.obj
 
 PLUGIN_DLL 	=	gdal_netCDF.dll
 

gdal/frmts/netcdf/netcdfdataset.cpp

@@ -51,6 +51,7 @@
 // Must be included after standard includes, otherwise VS2015 fails when
 // including <ctime>
 #include "netcdfdataset.h"
+#include "netcdfsg.h"
 #include "netcdfuffd.h"
 
 #include "cpl_conv.h"
@@ -2103,9 +2104,10 @@ netCDFDataset::netCDFDataset() :
     eFormat(NCDF_FORMAT_NONE),
     bIsGdalFile(false),
     bIsGdalCfFile(false),
-
     pszCFProjection(nullptr),
     pszCFCoordinates(nullptr),
+    nCFVersion(1.6),
+    bSGSupport(false),
     eMultipleLayerBehaviour(SINGLE_LAYER),
 
     // projection/GT.
@@ -2134,6 +2136,7 @@ netCDFDataset::netCDFDataset() :
     bSignedData(true),
     nLayers(0),
     papoLayers(nullptr)
+
 {
     // Projection/GT.
     adfGeoTransform[0] = 0.0;
@@ -2413,7 +2416,7 @@ static bool IsDifferenceBelow(double dfA, double dfB, double dfError)
 /*                      SetProjectionFromVar()                          */
 /************************************************************************/
 void netCDFDataset::SetProjectionFromVar( int nGroupId, int nVarId,
-                                          bool bReadSRSOnly )
+                                          bool bReadSRSOnly, const char * pszGivenGM)
 {
     bool bGotGeogCS = false;
     bool bGotCfSRS = false;
@@ -2456,7 +2459,8 @@ void netCDFDataset::SetProjectionFromVar( int nGroupId, int nVarId,
     }
 
     // Look for grid_mapping metadata.
-    const char *pszValue = FetchAttr(nGroupId, nVarId, CF_GRD_MAPPING);
+    const char *pszValue = pszGivenGM != nullptr ? pszGivenGM
+                                                 : FetchAttr(nGroupId, nVarId, CF_GRD_MAPPING);
     char *pszGridMappingValue = CPLStrdup(pszValue ? pszValue : "");
 
     if( !EQUAL(pszGridMappingValue, "") )
@@ -3872,6 +3876,12 @@ void netCDFDataset::SetProjectionFromVar( int nGroupId, int nVarId,
 #endif
 }
 
+void netCDFDataset::SetProjectionFromVar( int nGroupId, int nVarId,
+                                          bool bReadSRSOnly )
+{
+    SetProjectionFromVar(nGroupId, nVarId, bReadSRSOnly, nullptr);
+}
+
 int netCDFDataset::ProcessCFGeolocation( int nGroupId, int nVarId )
 {
     bool bAddGeoloc = false;
@@ -7127,6 +7137,19 @@ GDALDataset *netCDFDataset::Open( GDALOpenInfo *poOpenInfo )
         }
     }
 
+    // Figure out whether or not the listed dataset has support for simple geometries (CF-1.8)
+    poDS->nCFVersion = nccfdriver::getCFVersion(cdfid);
+    if(poDS->nCFVersion >= 1.8)
+    {
+        poDS->bSGSupport = true;
+    }
+    else poDS->bSGSupport = false;
+
+    if(poDS->bSGSupport)
+    {
+         poDS->DetectAndFillSGLayers(cdfid);
+    }
+
     char szConventions[NC_MAX_NAME + 1];
     szConventions[0] = '\0';
     nc_type nAttype = NC_NAT;
@@ -10909,7 +10932,7 @@ CPLErr netCDFDataset::FilterVars( int nCdfId, bool bKeepRasters,
         *pnRasterVars += nRasterVars;
     }
 
-    if( !anPotentialVectorVarID.empty() && bKeepVectors )
+    if( !anPotentialVectorVarID.empty() && bKeepVectors && !bSGSupport)
     {
         // Take the dimension that is referenced the most times.
         if( !(oMapDimIdToCount.size() == 1 ||

gdal/frmts/netcdf/netcdfdataset.h

@@ -39,6 +39,7 @@
 #include "gdal_pam.h"
 #include "gdal_priv.h"
 #include "netcdf.h"
+#include "netcdfsg.h"
 #include "ogr_spatialref.h"
 #include "ogrsf_frmts.h"
 #include "netcdfuffd.h"
@@ -49,6 +50,11 @@
 #define ENABLE_NCDUMP
 #endif
 
+#if CPL_IS_LSB
+#define PLATFORM_HEADER 1
+#else
+#define PLATFORM_HEADER 0
+#endif
 
 /************************************************************************/
 /* ==================================================================== */
@@ -84,6 +90,7 @@
 
 /* NETCDF driver defs */
 static const size_t NCDF_MAX_STR_LEN = 8192;
+#define NCDF_CONVENTIONS     "Conventions"
 #define NCDF_CONVENTIONS_CF_V1_5  "CF-1.5"
 #define NCDF_CONVENTIONS_CF_V1_6  "CF-1.6"
 #define NCDF_SPATIAL_REF     "spatial_ref"
@@ -250,6 +257,19 @@ static const int NCDF_DEFLATE_LEVEL    = 1;  /* best time/size ratio */
 #define CF_PP_GRID_NORTH_POLE_LATITUDE  "grid_north_pole_latitude"
 #define CF_PP_NORTH_POLE_GRID_LONGITUDE "north_pole_grid_longitude"
 
+/* Simple Geometries Special Names from CF-1.8 Draft - Chapter 7 section Geometries */
+#define CF_SG_GEOMETRY               "geometry"
+#define CF_SG_GEOMETRY_DIMENSION     "geometry_dimension"
+#define CF_SG_GEOMETRY_TYPE          "geometry_type"
+#define CF_SG_INTERIOR_RING          "interior_ring"
+#define CF_SG_NODES                  "nodes"
+#define CF_SG_NODE_COORDINATES       "node_coordinates"
+#define CF_SG_NODE_COUNT             "node_count"
+#define CF_SG_PART_NODE_COUNT        "part_node_count"
+#define CF_SG_TYPE_LINE              "line"
+#define CF_SG_TYPE_POINT             "point"
+#define CF_SG_TYPE_POLY              "polygon"
+
 /* -------------------------------------------------------------------- */
 /*         CF-1 Coordinate Type Naming (Chapter 4.  Coordinate Types )  */
 /* -------------------------------------------------------------------- */
@@ -808,6 +828,8 @@ class netCDFDataset final: public GDALPamDataset
     bool          bIsGdalCfFile; /* was this file created by the (new) CF-compliant driver? */
     char         *pszCFProjection;
     const char   *pszCFCoordinates;
+    double        nCFVersion;
+    bool          bSGSupport;
     MultipleLayerBehaviour eMultipleLayerBehaviour;
     std::vector<netCDFDataset*> apoVectorDatasets;
 
@@ -868,6 +890,7 @@ class netCDFDataset final: public GDALPamDataset
 
     void  CreateSubDatasetList( int nGroupId );
 
+    void  SetProjectionFromVar( int nGroupId, int nVarId, bool bReadSRSOnly, const char * pszGivenGM);
     void  SetProjectionFromVar( int nGroupId, int nVarId, bool bReadSRSOnly );
 
     int ProcessCFGeolocation( int nGroupId, int nVarId );
@@ -889,6 +912,10 @@ class netCDFDataset final: public GDALPamDataset
                                   int nVarXId, int nVarYId, int nVarZId,
                                   int nProfileDimId, int nParentIndexVarID,
                                   bool bKeepRasters );
+
+    CPLErr DetectAndFillSGLayers( int ncid );
+    CPLErr LoadSGVarIntoLayer( int ncid, int nc_basevarId );
+
   protected:
 
     CPLXMLNode *SerializeToXML( const char *pszVRTPath ) override;
@@ -944,6 +971,7 @@ class netCDFDataset final: public GDALPamDataset
 
 class netCDFLayer final: public OGRLayer
 {
+        friend class netCDFDataset;
         typedef union
         {
             signed char chVal;
@@ -994,6 +1022,10 @@ class netCDFLayer final: public OGRLayer
         nc_type         m_nWKTNCDFType;
         CPLString       m_osCoordinatesValue;
         std::vector<FieldDesc> m_aoFieldDesc;
+        std::vector<std::unique_ptr<OGRFeature>> m_sgFeatureList;
+        std::vector<std::unique_ptr<OGRFeature>>::iterator m_sgFeatItr;
+        bool            m_sgItrInit;
+        bool            m_HasCFSG1_8;
         int             m_nCurFeatureId;
         CPLString       m_osGridMapping;
         bool            m_bWriteGDALTags;
@@ -1019,9 +1051,10 @@ class netCDFLayer final: public OGRLayer
         void            GetNoDataValueForFloat( int nVarId, NCDFNoDataUnion* puNoData );
         void            GetNoDataValueForDouble( int nVarId, NCDFNoDataUnion* puNoData );
         void            GetNoDataValue( int nVarId, nc_type nVarType, NCDFNoDataUnion* puNoData );
-        bool            FillFeatureFromVar(OGRFeature* poFeature, int nMainDimId, size_t nIndex);
         bool            FillVarFromFeature(OGRFeature* poFeature, int nMainDimId, size_t nIndex);
 
+    protected:
+        bool            FillFeatureFromVar(OGRFeature* poFeature, int nMainDimId, size_t nIndex);
     public:
                 netCDFLayer(netCDFDataset* poDS,
                             int nLayerCDFId,
@@ -1036,6 +1069,8 @@ class netCDFLayer final: public OGRLayer
         void            SetWKTGeometryField(const char* pszWKTVarName);
         void            SetGridMapping(const char* pszGridMapping);
         void            SetProfile(int nProfileDimID, int nParentIndexVarID);
+        void            AddSimpleGeometryFeature(OGRFeature* sg) { this->m_sgFeatureList.push_back(std::unique_ptr<OGRFeature>(sg)); }
+        void            EnableSGBypass() { this-> m_HasCFSG1_8 = true; }
         bool            AddField(int nVarId);
 
         int             GetCDFID() const { return m_nLayerCDFId; }