Additions to netCDF Driver to allow Reading from CF-1.8 Encoded Geometries (#1287)

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),
+    nCFMinorVersion(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;
@@ -7127,6 +7130,19 @@ GDALDataset *netCDFDataset::Open( GDALOpenInfo *poOpenInfo )
         }
     }
 
+    // Figure out whether or not the listed dataset has support for simple geometries (CF-1.8)
+    poDS->nCFMinorVersion = nccfdriver::getCFMinorVersion(cdfid);
+    if(poDS->nCFMinorVersion >= 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;
@@ -10916,6 +10932,7 @@ CPLErr netCDFDataset::FilterVars( int nCdfId, bool bKeepRasters,
             (EQUAL(osFeatureType, "profile") &&
              oMapDimIdToCount.size() == 2 && nProfileDimId >= 0)) )
         {
+	if(!bSGSupport)
             CPLError(CE_Warning, CPLE_AppDefined,
                     "The dataset has several variables that could be "
                     "identified as vector fields, but not all share the same "
@@ -10954,7 +10971,7 @@ CPLErr netCDFDataset::CreateGrpVectorLayers( int nCdfId,
                                              int nVarXId, int nVarYId, int nVarZId,
                                              int nProfileDimId,
                                              int nParentIndexVarID,
-                                             bool bKeepRasters )
+					      bool bKeepRasters )
 {
     char *pszGroupName = nullptr;
     NCDFGetGroupFullName(nCdfId, &pszGroupName);

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"
@@ -84,6 +85,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 +252,18 @@ 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_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 +822,8 @@ class netCDFDataset final: public GDALPamDataset
     bool          bIsGdalCfFile; /* was this file created by the (new) CF-compliant driver? */
     char         *pszCFProjection;
     const char   *pszCFCoordinates;
+    int          nCFMinorVersion;
+    bool         bSGSupport;
     MultipleLayerBehaviour eMultipleLayerBehaviour;
     std::vector<netCDFDataset*> apoVectorDatasets;
 
@@ -889,6 +905,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 +964,7 @@ class netCDFDataset final: public GDALPamDataset
 
 class netCDFLayer final: public OGRLayer
 {
+	friend class netCDFDataset;
         typedef union
         {
             signed char chVal;
@@ -994,6 +1015,10 @@ class netCDFLayer final: public OGRLayer
         nc_type         m_nWKTNCDFType;
         CPLString       m_osCoordinatesValue;
         std::vector<FieldDesc> m_aoFieldDesc;
+	std::vector<OGRFeature*> m_sgFeatureList;
+	std::vector<OGRFeature*>::iterator m_sgFeatItr;
+	bool		m_sgItrInit;
+	bool		m_HasCFSG1_8;
         int             m_nCurFeatureId;
         CPLString       m_osGridMapping;
         bool            m_bWriteGDALTags;
@@ -1019,9 +1044,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 +1062,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(sg); }
+	void		EnableSGBypass() { this-> m_HasCFSG1_8 = true; }
         bool            AddField(int nVarId);
 
         int             GetCDFID() const { return m_nLayerCDFId; }

gdal/frmts/netcdf/netcdflayer.cpp

@@ -27,6 +27,7 @@
  ****************************************************************************/
 
 #include "netcdfdataset.h"
+#include "netcdfsg.h"
 #include "cpl_time.h"
 
 CPL_CVSID("$Id$")
@@ -59,6 +60,8 @@ netCDFLayer::netCDFLayer(netCDFDataset *poDS,
         m_nWKTMaxWidthDimId(-1),
         m_nWKTVarID(-1),
         m_nWKTNCDFType(NC_NAT),
+	m_sgItrInit(false),
+	m_HasCFSG1_8(false),
         m_nCurFeatureId(1),
         m_bWriteGDALTags(true),
         m_bUseStringInNC4(true),
@@ -83,7 +86,10 @@ netCDFLayer::netCDFLayer(netCDFDataset *poDS,
 /*                           ~netCDFLayer()                             */
 /************************************************************************/
 
-netCDFLayer::~netCDFLayer() { m_poFeatureDefn->Release(); }
+netCDFLayer::~netCDFLayer()
+{ 
+	m_poFeatureDefn->Release();
+}
 
 /************************************************************************/
 /*                   netCDFWriteAttributesFromConf()                    */
@@ -163,6 +169,7 @@ bool netCDFLayer::Create(char **papszOptions,
                          "geometry type.");
             }
         }
+
         else if( EQUAL(pszFeatureType, "PROFILE") )
         {
             if( wkbFlatten(m_poFeatureDefn->GetGeomType()) != wkbPoint )
@@ -1062,6 +1069,10 @@ bool netCDFLayer::FillFeatureFromVar(OGRFeature *poFeature, int nMainDimId,
         }
     }
 
+    // For CF-1.8 simple geometry specifically 
+    // Only need fields to be set here
+    if( this-> m_HasCFSG1_8 ) return true;
+
     if( m_nXVarID >= 0 && m_nYVarID >= 0 &&
         (m_osProfileDimName.empty() || nMainDimId == m_nProfileDimID) )
     {
@@ -1173,8 +1184,20 @@ bool netCDFLayer::FillFeatureFromVar(OGRFeature *poFeature, int nMainDimId,
 
 OGRFeature *netCDFLayer::GetNextFeature()
 {
+    if(!this->m_sgItrInit && m_sgFeatureList.size() > 0)
+    {
+	this->m_sgFeatItr = m_sgFeatureList.begin();
+	this->m_sgItrInit = true;
+    }
+
+    if(this->m_sgItrInit && this->m_sgFeatItr != m_sgFeatureList.end())
+    {
+	return *(m_sgFeatItr++);
+    }
+
     while( true )
     {
+
         OGRFeature *poFeature = GetNextRawFeature();
         if( poFeature == nullptr )
             return nullptr;
@@ -2059,7 +2082,7 @@ OGRErr netCDFLayer::CreateField(OGRFieldDefn *poFieldDefn, int /* bApproxOK */)
             osVarName = CPLSPrintf("%s%d", poFieldDefn->GetNameRef(), i);
             status = nc_inq_varid(m_nLayerCDFId, osVarName, &nVarID);
             if( status != NC_NOERR )
-                break;
+            break;
         }
         CPLDebug("netCDF", "Field %s is written in variable %s",
                  poFieldDefn->GetNameRef(), osVarName.c_str());