/****************************************************************************** * * Project: netCDF read/write Driver * Purpose: GDAL bindings over netCDF library. * Author: Even Rouault * ****************************************************************************** * Copyright (c) 2016, Even Rouault * * Permission is hereby granted, free of charge, to any person obtaining a * copy of this software and associated documentation files (the "Software"), * to deal in the Software without restriction, including without limitation * the rights to use, copy, modify, merge, publish, distribute, sublicense, * and/or sell copies of the Software, and to permit persons to whom the * Software is furnished to do so, subject to the following conditions: * * The above copyright notice and this permission notice shall be included * in all copies or substantial portions of the Software. * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER * DEALINGS IN THE SOFTWARE. ****************************************************************************/ #include "netcdfdataset.h" #include "cpl_time.h" CPL_CVSID("$Id: netcdflayer.cpp 7d078e0357d2998edfa713422e607cbadf77f9ff 2018-04-08 22:11:28 +0200 Even Rouault $") /************************************************************************/ /* netCDFLayer() */ /************************************************************************/ netCDFLayer::netCDFLayer(netCDFDataset *poDS, int nLayerCDFId, const char *pszName, OGRwkbGeometryType eGeomType, OGRSpatialReference *poSRS) : m_poDS(poDS), m_nLayerCDFId(nLayerCDFId), m_poFeatureDefn(new OGRFeatureDefn(pszName)), m_osRecordDimName("record"), m_nRecordDimID(-1), m_nDefaultWidth(10), m_bAutoGrowStrings(true), m_nDefaultMaxWidthDimId(-1), m_nXVarID(-1), m_nYVarID(-1), m_nZVarID(-1), m_nXVarNCDFType(NC_NAT), m_nYVarNCDFType(NC_NAT), m_nZVarNCDFType(NC_NAT), m_osWKTVarName("ogc_wkt"), m_nWKTMaxWidth(10000), m_nWKTMaxWidthDimId(-1), m_nWKTVarID(-1), m_nWKTNCDFType(NC_NAT), m_nCurFeatureId(1), m_bWriteGDALTags(true), m_bUseStringInNC4(true), m_bNCDumpCompat(true), m_nProfileDimID(-1), m_nProfileVarID(-1), m_bProfileVarUnlimited(false), m_nParentIndexVarID(-1), m_poLayerConfig(nullptr) { m_uXVarNoData.nVal64 = 0; m_uYVarNoData.nVal64 = 0; m_uZVarNoData.nVal64 = 0; m_poFeatureDefn->SetGeomType(eGeomType); if( eGeomType != wkbNone ) m_poFeatureDefn->GetGeomFieldDefn(0)->SetSpatialRef(poSRS); m_poFeatureDefn->Reference(); SetDescription(pszName); } /************************************************************************/ /* ~netCDFLayer() */ /************************************************************************/ netCDFLayer::~netCDFLayer() { m_poFeatureDefn->Release(); } /************************************************************************/ /* netCDFWriteAttributesFromConf() */ /************************************************************************/ static void netCDFWriteAttributesFromConf( int cdfid, int varid, const std::vector &aoAttributes) { for(size_t i = 0; i < aoAttributes.size(); i++) { const netCDFWriterConfigAttribute &oAtt = aoAttributes[i]; int status = NC_NOERR; if( oAtt.m_osValue.empty() ) { int attid = -1; status = nc_inq_attid(cdfid, varid, oAtt.m_osName, &attid); if( status == NC_NOERR ) status = nc_del_att(cdfid, varid, oAtt.m_osName); else status = NC_NOERR; } else if( EQUAL(oAtt.m_osType, "string") ) { status = nc_put_att_text(cdfid, varid, oAtt.m_osName, oAtt.m_osValue.size(), oAtt.m_osValue); } else if( EQUAL(oAtt.m_osType, "integer") ) { int nVal = atoi(oAtt.m_osValue); status = nc_put_att_int(cdfid, varid, oAtt.m_osName, NC_INT, 1, &nVal); } else if( EQUAL(oAtt.m_osType, "double") ) { double dfVal = CPLAtof(oAtt.m_osValue); status = nc_put_att_double(cdfid, varid, oAtt.m_osName, NC_DOUBLE, 1, &dfVal); } NCDF_ERR(status); } } /************************************************************************/ /* Create() */ /************************************************************************/ bool netCDFLayer::Create(char **papszOptions, const netCDFWriterConfigLayer *poLayerConfig) { m_osRecordDimName = CSLFetchNameValueDef(papszOptions, "RECORD_DIM_NAME", m_osRecordDimName.c_str()); m_bAutoGrowStrings = CPL_TO_BOOL(CSLFetchBoolean(papszOptions, "AUTOGROW_STRINGS", TRUE)); m_nDefaultWidth = atoi( CSLFetchNameValueDef(papszOptions, "STRING_DEFAULT_WIDTH", CPLSPrintf("%d", m_bAutoGrowStrings ? 10 : 80))); m_bWriteGDALTags = CPL_TO_BOOL( CSLFetchBoolean(m_poDS->papszCreationOptions, "WRITE_GDAL_TAGS", TRUE)); m_bUseStringInNC4 = CPL_TO_BOOL(CSLFetchBoolean(papszOptions, "USE_STRING_IN_NC4", TRUE)); m_bNCDumpCompat = CPL_TO_BOOL(CSLFetchBoolean(papszOptions, "NCDUMP_COMPAT", TRUE)); std::vector > aoAutoVariables; const char *pszFeatureType = CSLFetchNameValue(papszOptions, "FEATURE_TYPE"); if( pszFeatureType != nullptr ) { if( EQUAL(pszFeatureType, "POINT") ) { if( wkbFlatten(m_poFeatureDefn->GetGeomType()) != wkbPoint ) { CPLError(CE_Warning, CPLE_NotSupported, "FEATURE_TYPE=POINT only supported for Point layer " "geometry type."); } } else if( EQUAL(pszFeatureType, "PROFILE") ) { if( wkbFlatten(m_poFeatureDefn->GetGeomType()) != wkbPoint ) { CPLError(CE_Warning, CPLE_NotSupported, "FEATURE_TYPE=PROFILE only supported for Point layer " "geometry type."); } else { // Cf http://cfconventions.org/cf-conventions/v1.6.0/cf-conventions.html#_indexed_ragged_array_representation_of_profiles m_osProfileDimName = CSLFetchNameValueDef( papszOptions, "PROFILE_DIM_NAME", "profile"); m_osProfileVariables = CSLFetchNameValueDef(papszOptions, "PROFILE_VARIABLES", ""); const char *pszProfileInitSize = CSLFetchNameValueDef( papszOptions, "PROFILE_DIM_INIT_SIZE", (m_poDS->eFormat == NCDF_FORMAT_NC4) ? "UNLIMITED" : "100"); m_bProfileVarUnlimited = EQUAL(pszProfileInitSize, "UNLIMITED"); size_t nProfileInitSize = m_bProfileVarUnlimited ? NC_UNLIMITED : static_cast(atoi(pszProfileInitSize)); int status = nc_def_dim(m_nLayerCDFId, m_osProfileDimName, nProfileInitSize, &m_nProfileDimID); NCDF_ERR(status); if( status != NC_NOERR ) return false; status = nc_def_var(m_nLayerCDFId, m_osProfileDimName, NC_INT, 1, &m_nProfileDimID, &m_nProfileVarID); NCDF_ERR(status); if( status != NC_NOERR ) return false; aoAutoVariables.push_back(std::pair( m_osProfileDimName, m_nProfileVarID)); status = nc_put_att_text(m_nLayerCDFId, m_nProfileVarID, "cf_role", strlen("profile_id"), "profile_id"); NCDF_ERR(status); } } else if( !EQUAL(pszFeatureType, "AUTO") ) { CPLError(CE_Warning, CPLE_NotSupported, "FEATURE_TYPE=%s not supported.", pszFeatureType); } } int status; if( m_bWriteGDALTags ) { status = nc_put_att_text(m_nLayerCDFId, NC_GLOBAL, "ogr_layer_name", strlen(m_poFeatureDefn->GetName()), m_poFeatureDefn->GetName()); NCDF_ERR(status); } status = nc_def_dim(m_nLayerCDFId, m_osRecordDimName, NC_UNLIMITED, &m_nRecordDimID); NCDF_ERR(status); if( status != NC_NOERR ) return false; if( !m_osProfileDimName.empty() ) { status = nc_def_var(m_nLayerCDFId, "parentIndex", NC_INT, 1, &m_nRecordDimID, &m_nParentIndexVarID); NCDF_ERR(status); if( status != NC_NOERR ) return false; aoAutoVariables.push_back( std::pair("parentIndex", m_nParentIndexVarID)); status = nc_put_att_text(m_nLayerCDFId, m_nParentIndexVarID, CF_LNG_NAME, strlen("index of profile"), "index of profile"); NCDF_ERR(status); status = nc_put_att_text( m_nLayerCDFId, m_nParentIndexVarID, "instance_dimension", m_osProfileDimName.size(), m_osProfileDimName.c_str()); NCDF_ERR(status); } OGRSpatialReference *poSRS = nullptr; if( m_poFeatureDefn->GetGeomFieldCount() ) poSRS = m_poFeatureDefn->GetGeomFieldDefn(0)->GetSpatialRef(); if( wkbFlatten(m_poFeatureDefn->GetGeomType()) == wkbPoint ) { const int nPointDim = !m_osProfileDimName.empty() ? m_nProfileDimID : m_nRecordDimID; const bool bIsGeographic = (poSRS == nullptr || poSRS->IsGeographic()); const char *pszXVarName = bIsGeographic ? CF_LONGITUDE_VAR_NAME : CF_PROJ_X_VAR_NAME; status = nc_def_var(m_nLayerCDFId, pszXVarName, NC_DOUBLE, 1, &nPointDim, &m_nXVarID); NCDF_ERR(status); if( status != NC_NOERR ) { return false; } const char *pszYVarName = bIsGeographic ? CF_LATITUDE_VAR_NAME : CF_PROJ_Y_VAR_NAME; status = nc_def_var(m_nLayerCDFId, pszYVarName, NC_DOUBLE, 1, &nPointDim, &m_nYVarID); NCDF_ERR(status); if( status != NC_NOERR ) { return false; } aoAutoVariables.push_back( std::pair(pszXVarName, m_nXVarID)); aoAutoVariables.push_back( std::pair(pszYVarName, m_nYVarID)); m_nXVarNCDFType = NC_DOUBLE; m_nYVarNCDFType = NC_DOUBLE; m_uXVarNoData.dfVal = NC_FILL_DOUBLE; m_uYVarNoData.dfVal = NC_FILL_DOUBLE; m_osCoordinatesValue = pszXVarName; m_osCoordinatesValue += " "; m_osCoordinatesValue += pszYVarName; if( poSRS == nullptr || poSRS->IsGeographic() ) { NCDFWriteLonLatVarsAttributes(m_nLayerCDFId, m_nXVarID, m_nYVarID); } else if (poSRS != nullptr && poSRS->IsProjected() ) { NCDFWriteXYVarsAttributes(m_nLayerCDFId, m_nXVarID, m_nYVarID, poSRS); } if( m_poFeatureDefn->GetGeomType() == wkbPoint25D ) { const char *pszZVarName = "z"; status = nc_def_var(m_nLayerCDFId, pszZVarName, NC_DOUBLE, 1, &m_nRecordDimID, &m_nZVarID); NCDF_ERR(status); if( status != NC_NOERR ) { return false; } aoAutoVariables.push_back( std::pair(pszZVarName, m_nZVarID)); m_nZVarNCDFType = NC_DOUBLE; m_uZVarNoData.dfVal = NC_FILL_DOUBLE; status = nc_put_att_text(m_nLayerCDFId, m_nZVarID, CF_LNG_NAME, strlen("z coordinate"), "z coordinate"); NCDF_ERR(status); status = nc_put_att_text(m_nLayerCDFId, m_nZVarID, CF_STD_NAME, strlen("height"), "height"); NCDF_ERR(status); status = nc_put_att_text(m_nLayerCDFId, m_nZVarID, CF_AXIS, strlen("Z"), "Z"); NCDF_ERR(status); status = nc_put_att_text(m_nLayerCDFId, m_nZVarID, CF_UNITS, strlen("m"), "m"); NCDF_ERR(status); m_osCoordinatesValue += " "; m_osCoordinatesValue += pszZVarName; } const char *pszFeatureTypeVal = !m_osProfileDimName.empty() ? "profile" : "point"; status = nc_put_att_text(m_nLayerCDFId, NC_GLOBAL, "featureType", strlen(pszFeatureTypeVal), pszFeatureTypeVal); NCDF_ERR(status); } else if( m_poFeatureDefn->GetGeomType() != wkbNone ) { #ifdef NETCDF_HAS_NC4 if( m_poDS->eFormat == NCDF_FORMAT_NC4 && m_bUseStringInNC4 ) { m_nWKTNCDFType = NC_STRING; status = nc_def_var(m_nLayerCDFId, m_osWKTVarName.c_str(), NC_STRING, 1, &m_nRecordDimID, &m_nWKTVarID); } else #endif { m_nWKTNCDFType = NC_CHAR; m_nWKTMaxWidth = atoi(CSLFetchNameValueDef( papszOptions, "WKT_DEFAULT_WIDTH", CPLSPrintf("%d", m_bAutoGrowStrings ? 1000 : 10000))); status = nc_def_dim(m_nLayerCDFId, CPLSPrintf("%s_max_width", m_osWKTVarName.c_str()), m_nWKTMaxWidth, &m_nWKTMaxWidthDimId); NCDF_ERR(status); if( status != NC_NOERR ) { return false; } int anDims[2] = { m_nRecordDimID, m_nWKTMaxWidthDimId }; status = nc_def_var(m_nLayerCDFId, m_osWKTVarName.c_str(), NC_CHAR, 2, anDims, &m_nWKTVarID); } NCDF_ERR(status); if( status != NC_NOERR ) { return false; } aoAutoVariables.push_back( std::pair(m_osWKTVarName, m_nWKTVarID)); status = nc_put_att_text(m_nLayerCDFId, m_nWKTVarID, CF_LNG_NAME, strlen("Geometry as ISO WKT"), "Geometry as ISO WKT"); NCDF_ERR(status); // nc_put_att_text(m_nLayerCDFId, m_nWKTVarID, CF_UNITS, // strlen("none"), "none"); if( m_bWriteGDALTags ) { status = nc_put_att_text(m_nLayerCDFId, NC_GLOBAL, "ogr_geometry_field", m_osWKTVarName.size(), m_osWKTVarName.c_str()); NCDF_ERR(status); CPLString osGeometryType = OGRToOGCGeomType(m_poFeatureDefn->GetGeomType()); if( wkbHasZ(m_poFeatureDefn->GetGeomType()) ) osGeometryType += " Z"; status = nc_put_att_text(m_nLayerCDFId, NC_GLOBAL, "ogr_layer_type", osGeometryType.size(), osGeometryType.c_str()); NCDF_ERR(status); } } if( poSRS != nullptr ) { char *pszCFProjection = nullptr; int nSRSVarId = NCDFWriteSRSVariable( m_nLayerCDFId, poSRS, &pszCFProjection, m_bWriteGDALTags); if( nSRSVarId < 0 ) return false; if( pszCFProjection != nullptr ) { aoAutoVariables.push_back( std::pair(pszCFProjection, nSRSVarId)); m_osGridMapping = pszCFProjection; CPLFree(pszCFProjection); } if( m_nWKTVarID >= 0 && !m_osGridMapping.empty() ) { status = nc_put_att_text(m_nLayerCDFId, m_nWKTVarID, CF_GRD_MAPPING, m_osGridMapping.size(), m_osGridMapping.c_str()); NCDF_ERR(status); } } if( m_poDS->oWriterConfig.m_bIsValid ) { m_poLayerConfig = poLayerConfig; netCDFWriteAttributesFromConf(m_nLayerCDFId, NC_GLOBAL, m_poDS->oWriterConfig.m_aoAttributes); if( poLayerConfig != nullptr ) { netCDFWriteAttributesFromConf(m_nLayerCDFId, NC_GLOBAL, poLayerConfig->m_aoAttributes); } for( size_t i = 0; i < aoAutoVariables.size(); i++ ) { const netCDFWriterConfigField *poConfig = nullptr; CPLString osLookup = "__" + aoAutoVariables[i].first; std::map::const_iterator oIter; if( m_poLayerConfig != nullptr && (oIter = m_poLayerConfig->m_oFields.find(osLookup)) != m_poLayerConfig->m_oFields.end() ) { poConfig = &(oIter->second); } else if( (oIter = m_poDS->oWriterConfig.m_oFields.find(osLookup)) != m_poDS->oWriterConfig.m_oFields.end() ) { poConfig = &(oIter->second); } if( poConfig != nullptr ) { netCDFWriteAttributesFromConf(m_nLayerCDFId, aoAutoVariables[i].second, poConfig->m_aoAttributes); } } } return true; } /************************************************************************/ /* SetRecordDimID() */ /************************************************************************/ void netCDFLayer::SetRecordDimID(int nRecordDimID) { m_nRecordDimID = nRecordDimID; char szTemp[NC_MAX_NAME + 1]; szTemp[0] = 0; int status = nc_inq_dimname(m_nLayerCDFId, m_nRecordDimID, szTemp); NCDF_ERR(status); m_osRecordDimName = szTemp; } /************************************************************************/ /* GetFillValue() */ /************************************************************************/ CPLErr netCDFLayer::GetFillValue( int nVarId, char **ppszValue ) { if( NCDFGetAttr(m_nLayerCDFId, nVarId, _FillValue, ppszValue) == CE_None ) return CE_None; return NCDFGetAttr(m_nLayerCDFId, nVarId, "missing_value", ppszValue); } CPLErr netCDFLayer::GetFillValue( int nVarId, double *pdfValue ) { if( NCDFGetAttr(m_nLayerCDFId, nVarId, _FillValue, pdfValue) == CE_None ) return CE_None; return NCDFGetAttr(m_nLayerCDFId, nVarId, "missing_value", pdfValue); } /************************************************************************/ /* GetNoDataValueForFloat() */ /************************************************************************/ void netCDFLayer::GetNoDataValueForFloat( int nVarId, NCDFNoDataUnion *puNoData ) { double dfValue; if( GetFillValue(nVarId, &dfValue) == CE_None ) puNoData->fVal = static_cast(dfValue); else puNoData->fVal = NC_FILL_FLOAT; } /************************************************************************/ /* GetNoDataValueForDouble() */ /************************************************************************/ void netCDFLayer::GetNoDataValueForDouble( int nVarId, NCDFNoDataUnion *puNoData ) { double dfValue; if( GetFillValue(nVarId, &dfValue) == CE_None ) puNoData->dfVal = dfValue; else puNoData->dfVal = NC_FILL_DOUBLE; } /************************************************************************/ /* GetNoDataValue() */ /************************************************************************/ void netCDFLayer::GetNoDataValue( int nVarId, nc_type nVarType, NCDFNoDataUnion *puNoData ) { if( nVarType == NC_DOUBLE ) GetNoDataValueForDouble(nVarId, puNoData); else if( nVarType == NC_FLOAT ) GetNoDataValueForFloat(nVarId, puNoData); } /************************************************************************/ /* SetXYZVars() */ /************************************************************************/ void netCDFLayer::SetXYZVars(int nXVarId, int nYVarId, int nZVarId) { m_nXVarID = nXVarId; m_nYVarID = nYVarId; m_nZVarID = nZVarId; nc_inq_vartype(m_nLayerCDFId, m_nXVarID, &m_nXVarNCDFType); nc_inq_vartype(m_nLayerCDFId, m_nYVarID, &m_nYVarNCDFType); if( (m_nXVarNCDFType != NC_FLOAT && m_nXVarNCDFType != NC_DOUBLE) || (m_nYVarNCDFType != NC_FLOAT && m_nYVarNCDFType != NC_DOUBLE) ) { CPLError(CE_Warning, CPLE_NotSupported, "X or Y variable of type X=%d,Y=%d not handled", m_nXVarNCDFType, m_nYVarNCDFType); m_nXVarID = -1; m_nYVarID = -1; } if( m_nZVarID >= 0 ) { nc_inq_vartype(m_nLayerCDFId, m_nZVarID, &m_nZVarNCDFType); if( m_nZVarNCDFType != NC_FLOAT && m_nZVarNCDFType != NC_DOUBLE ) { CPLError(CE_Warning, CPLE_NotSupported, "Z variable of type %d not handled", m_nZVarNCDFType); m_nZVarID = -1; } } if( m_nXVarID >= 0 && m_nYVarID >= 0 ) { char szVarName[NC_MAX_NAME + 1]; szVarName[0] = '\0'; CPL_IGNORE_RET_VAL(nc_inq_varname(m_nLayerCDFId, m_nXVarID, szVarName)); m_osCoordinatesValue = szVarName; szVarName[0] = '\0'; CPL_IGNORE_RET_VAL(nc_inq_varname(m_nLayerCDFId, m_nYVarID, szVarName)); m_osCoordinatesValue += " "; m_osCoordinatesValue += szVarName; if( m_nZVarID >= 0 ) { szVarName[0] = '\0'; CPL_IGNORE_RET_VAL( nc_inq_varname(m_nLayerCDFId, m_nZVarID, szVarName)); m_osCoordinatesValue += " "; m_osCoordinatesValue += szVarName; } } if( m_nXVarID >= 0 ) GetNoDataValue(m_nXVarID, m_nXVarNCDFType, &m_uXVarNoData); if( m_nYVarID >= 0 ) GetNoDataValue(m_nYVarID, m_nYVarNCDFType, &m_uYVarNoData); if( m_nZVarID >= 0 ) GetNoDataValue(m_nZVarID, m_nZVarNCDFType, &m_uZVarNoData); } /************************************************************************/ /* SetWKTGeometryField() */ /************************************************************************/ void netCDFLayer::SetWKTGeometryField(const char *pszWKTVarName) { m_nWKTVarID = -1; nc_inq_varid(m_nLayerCDFId, pszWKTVarName, &m_nWKTVarID); if( m_nWKTVarID < 0 ) return; int nd; nc_inq_varndims(m_nLayerCDFId, m_nWKTVarID, &nd); nc_inq_vartype(m_nLayerCDFId, m_nWKTVarID, &m_nWKTNCDFType); #ifdef NETCDF_HAS_NC4 if( nd == 1 && m_nWKTNCDFType == NC_STRING ) { int nDimID; if( nc_inq_vardimid(m_nLayerCDFId, m_nWKTVarID, &nDimID ) != NC_NOERR || nDimID != m_nRecordDimID ) { m_nWKTVarID = -1; return; } } else #endif if (nd == 2 && m_nWKTNCDFType == NC_CHAR ) { int anDimIds [2] = { -1, -1 }; size_t nLen = 0; if( nc_inq_vardimid(m_nLayerCDFId, m_nWKTVarID, anDimIds) != NC_NOERR || anDimIds[0] != m_nRecordDimID || nc_inq_dimlen(m_nLayerCDFId, anDimIds[1], &nLen) != NC_NOERR ) { m_nWKTVarID = -1; return; } m_nWKTMaxWidth = static_cast(nLen); m_nWKTMaxWidthDimId = anDimIds[1]; } else { m_nWKTVarID = -1; return; } m_osWKTVarName = pszWKTVarName; } /************************************************************************/ /* SetGridMapping() */ /************************************************************************/ void netCDFLayer::SetGridMapping(const char *pszGridMapping) { m_osGridMapping = pszGridMapping; } /************************************************************************/ /* SetProfile() */ /************************************************************************/ void netCDFLayer::SetProfile(int nProfileDimID, int nParentIndexVarID) { m_nProfileDimID = nProfileDimID; m_nParentIndexVarID = nParentIndexVarID; if( m_nProfileDimID >= 0 ) { char szTemp[NC_MAX_NAME + 1]; szTemp[0] = 0; int status = nc_inq_dimname(m_nLayerCDFId, m_nProfileDimID, szTemp); NCDF_ERR(status); m_osProfileDimName = szTemp; nc_inq_varid(m_nLayerCDFId, m_osProfileDimName, &m_nProfileVarID); m_bProfileVarUnlimited = NCDFIsUnlimitedDim( m_poDS->eFormat == NCDF_FORMAT_NC4, m_nLayerCDFId, m_nProfileVarID); } } /************************************************************************/ /* ResetReading() */ /************************************************************************/ void netCDFLayer::ResetReading() { m_nCurFeatureId = 1; } /************************************************************************/ /* Get1DVarAsDouble() */ /************************************************************************/ double netCDFLayer::Get1DVarAsDouble( int nVarId, nc_type nVarType, size_t nIndex, NCDFNoDataUnion noDataVal, bool *pbIsNoData ) { double dfVal = 0; if( nVarType == NC_DOUBLE ) { nc_get_var1_double(m_nLayerCDFId, nVarId, &nIndex, &dfVal); if( pbIsNoData ) *pbIsNoData = dfVal == noDataVal.dfVal; } else if( nVarType == NC_FLOAT ) { float fVal = 0.f; nc_get_var1_float(m_nLayerCDFId, nVarId, &nIndex, &fVal); if( pbIsNoData ) *pbIsNoData = fVal == noDataVal.fVal; dfVal = fVal; } else if( pbIsNoData ) *pbIsNoData = true; return dfVal; } /************************************************************************/ /* GetNextRawFeature() */ /************************************************************************/ OGRFeature *netCDFLayer::GetNextRawFeature() { m_poDS->SetDefineMode(false); // In update mode, nc_get_varXXX() doesn't return error if we are // beyond the end of dimension size_t nDimLen = 0; nc_inq_dimlen(m_nLayerCDFId, m_nRecordDimID, &nDimLen); if( m_nCurFeatureId > static_cast(nDimLen) ) return nullptr; OGRFeature *poFeature = new OGRFeature(m_poFeatureDefn); if( m_nParentIndexVarID >= 0 ) { int nProfileIdx = 0; size_t nIdx = static_cast(m_nCurFeatureId - 1); int status = nc_get_var1_int(m_nLayerCDFId, m_nParentIndexVarID, &nIdx, &nProfileIdx); if( status == NC_NOERR && nProfileIdx >= 0 ) { nIdx = static_cast(nProfileIdx); FillFeatureFromVar(poFeature, m_nProfileDimID, nIdx); } } if( !FillFeatureFromVar(poFeature, m_nRecordDimID, static_cast(m_nCurFeatureId - 1)) ) { m_nCurFeatureId++; delete poFeature; return nullptr; } poFeature->SetFID(m_nCurFeatureId); m_nCurFeatureId++; return poFeature; } /************************************************************************/ /* FillFeatureFromVar() */ /************************************************************************/ bool netCDFLayer::FillFeatureFromVar(OGRFeature *poFeature, int nMainDimId, size_t nIndex) { size_t anIndex[2]; anIndex[0] = nIndex; anIndex[1] = 0; for( int i = 0; i < m_poFeatureDefn->GetFieldCount(); i++ ) { if( m_aoFieldDesc[i].nMainDimId != nMainDimId ) continue; switch( m_aoFieldDesc[i].nType ) { case NC_CHAR: { if( m_aoFieldDesc[i].nDimCount == 1 ) { char szVal[2] = { 0, 0 }; int status = nc_get_var1_text( m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, szVal); if( status != NC_NOERR ) { NCDF_ERR(status); continue; } poFeature->SetField(i, szVal); } else { size_t anCount[2]; anCount[0] = 1; anCount[1] = 0; nc_inq_dimlen(m_nLayerCDFId, m_aoFieldDesc[i].nSecDimId, &(anCount[1])); char *pszVal = (char *)CPLCalloc(1, anCount[1] + 1); int status = nc_get_vara_text(m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, anCount, pszVal); if( status != NC_NOERR ) { NCDF_ERR(status); CPLFree(pszVal); continue; } poFeature->SetField(i, pszVal); CPLFree(pszVal); } break; } #ifdef NETCDF_HAS_NC4 case NC_STRING: { char *pszVal = nullptr; int status = nc_get_var1_string( m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, &pszVal); if( status != NC_NOERR ) { NCDF_ERR(status); continue; } if( pszVal != nullptr ) { poFeature->SetField(i, pszVal); nc_free_string(1, &pszVal); } break; } #endif case NC_BYTE: { signed char chVal = 0; int status = nc_get_var1_schar( m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, &chVal); if( status != NC_NOERR ) { NCDF_ERR(status); continue; } if( chVal == m_aoFieldDesc[i].uNoData.chVal ) continue; poFeature->SetField(i, static_cast(chVal)); break; } #ifdef NETCDF_HAS_NC4 case NC_UBYTE: { unsigned char uchVal = 0; int status = nc_get_var1_uchar( m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, &uchVal); if( status != NC_NOERR ) { NCDF_ERR(status); continue; } if( uchVal == m_aoFieldDesc[i].uNoData.uchVal ) continue; poFeature->SetField(i, static_cast(uchVal)); break; } #endif case NC_SHORT: { short sVal = 0; int status = nc_get_var1_short( m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, &sVal); if( status != NC_NOERR ) { NCDF_ERR(status); continue; } if( sVal == m_aoFieldDesc[i].uNoData.sVal ) continue; poFeature->SetField(i, static_cast(sVal)); break; } #ifdef NETCDF_HAS_NC4 case NC_USHORT: { unsigned short usVal = 0; int status = nc_get_var1_ushort( m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, &usVal); if( status != NC_NOERR ) { NCDF_ERR(status); continue; } if( usVal == m_aoFieldDesc[i].uNoData.usVal ) continue; poFeature->SetField(i, static_cast(usVal)); break; } #endif case NC_INT: { int nVal = 0; int status = nc_get_var1_int(m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, &nVal); if( status != NC_NOERR ) { NCDF_ERR(status); continue; } if( nVal == m_aoFieldDesc[i].uNoData.nVal ) continue; if( m_poFeatureDefn->GetFieldDefn(i)->GetType() == OFTDate || m_poFeatureDefn->GetFieldDefn(i)->GetType() == OFTDateTime ) { struct tm brokendowntime; GIntBig nVal64 = static_cast(nVal); if( m_aoFieldDesc[i].bIsDays ) nVal64 *= 86400; CPLUnixTimeToYMDHMS(nVal64, &brokendowntime); poFeature->SetField( i, brokendowntime.tm_year + 1900, brokendowntime.tm_mon + 1, brokendowntime.tm_mday, brokendowntime.tm_hour, brokendowntime.tm_min, static_cast(brokendowntime.tm_sec), 0); } else { poFeature->SetField(i, nVal); } break; } #ifdef NETCDF_HAS_NC4 case NC_UINT: { unsigned int unVal = 0; // nc_get_var1_uint() doesn't work on old netCDF version when // the returned value is > INT_MAX // https://bugtracking.unidata.ucar.edu/browse/NCF-226 // nc_get_vara_uint() has not this bug size_t nCount = 1; int status = nc_get_vara_uint(m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, &nCount, &unVal); if( status != NC_NOERR ) { NCDF_ERR(status); continue; } if( unVal == m_aoFieldDesc[i].uNoData.unVal ) continue; poFeature->SetField(i, static_cast(unVal)); break; } #endif #ifdef NETCDF_HAS_NC4 case NC_INT64: { GIntBig nVal = 0; int status = nc_get_var1_longlong( m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, &nVal); if( status != NC_NOERR ) { NCDF_ERR(status); continue; } if( nVal == m_aoFieldDesc[i].uNoData.nVal64 ) continue; poFeature->SetField(i, nVal); break; } case NC_UINT64: { GUIntBig nVal = 0; int status = nc_get_var1_ulonglong( m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, &nVal); if( status != NC_NOERR ) { NCDF_ERR(status); continue; } if( nVal == m_aoFieldDesc[i].uNoData.unVal64 ) continue; poFeature->SetField(i, static_cast(nVal)); break; } #endif case NC_FLOAT: { float fVal = 0.f; int status = nc_get_var1_float( m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, &fVal); if( status != NC_NOERR ) { NCDF_ERR(status); continue; } if( fVal == m_aoFieldDesc[i].uNoData.fVal ) continue; poFeature->SetField(i, static_cast(fVal)); break; } case NC_DOUBLE: { double dfVal = 0.0; int status = nc_get_var1_double( m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, &dfVal); if( status != NC_NOERR ) { NCDF_ERR(status); continue; } if( dfVal == m_aoFieldDesc[i].uNoData.dfVal ) continue; if( m_poFeatureDefn->GetFieldDefn(i)->GetType() == OFTDate || m_poFeatureDefn->GetFieldDefn(i)->GetType() == OFTDateTime ) { if( m_aoFieldDesc[i].bIsDays ) dfVal *= 86400.0; struct tm brokendowntime; GIntBig nVal = static_cast(floor(dfVal)); CPLUnixTimeToYMDHMS(nVal, &brokendowntime); poFeature->SetField( i, brokendowntime.tm_year + 1900, brokendowntime.tm_mon + 1, brokendowntime.tm_mday, brokendowntime.tm_hour, brokendowntime.tm_min, static_cast(brokendowntime.tm_sec + (dfVal - nVal)), 0); } else { poFeature->SetField(i, dfVal); } break; } default: break; } } if( m_nXVarID >= 0 && m_nYVarID >= 0 && (m_osProfileDimName.empty() || nMainDimId == m_nProfileDimID) ) { bool bXIsNoData = false; const double dfX = Get1DVarAsDouble( m_nXVarID, m_nXVarNCDFType, anIndex[0], m_uXVarNoData, &bXIsNoData); bool bYIsNoData = false; const double dfY = Get1DVarAsDouble( m_nYVarID, m_nYVarNCDFType, anIndex[0], m_uYVarNoData, &bYIsNoData); if( !bXIsNoData && !bYIsNoData ) { OGRPoint *poPoint = nullptr; if( m_nZVarID >= 0 && m_osProfileDimName.empty() ) { bool bZIsNoData = false; const double dfZ = Get1DVarAsDouble(m_nZVarID, m_nZVarNCDFType, anIndex[0], m_uZVarNoData, &bZIsNoData); if( bZIsNoData ) poPoint = new OGRPoint(dfX, dfY); else poPoint = new OGRPoint(dfX, dfY, dfZ); } else poPoint = new OGRPoint(dfX, dfY); poPoint->assignSpatialReference(GetSpatialRef()); poFeature->SetGeometryDirectly(poPoint); } } else if( m_nXVarID >= 0 && m_nYVarID >= 0 && m_nZVarID >= 0 && !m_osProfileDimName.empty() && nMainDimId == m_nRecordDimID ) { OGRGeometry *poGeom = poFeature->GetGeometryRef(); if( poGeom != nullptr && wkbFlatten(poGeom->getGeometryType()) == wkbPoint ) { bool bZIsNoData = false; const double dfZ = Get1DVarAsDouble(m_nZVarID, m_nZVarNCDFType, anIndex[0], m_uZVarNoData, &bZIsNoData); if( !bZIsNoData ) poGeom->toPoint()->setZ(dfZ); } } else if( m_nWKTVarID >= 0 ) { char *pszWKT = nullptr; if( m_nWKTNCDFType == NC_CHAR ) { size_t anCount[2]; anCount[0] = 1; anCount[1] = m_nWKTMaxWidth; pszWKT = (char *)CPLCalloc(1, anCount[1] + 1); int status = nc_get_vara_text(m_nLayerCDFId, m_nWKTVarID, anIndex, anCount, pszWKT); if( status == NC_EINVALCOORDS || status == NC_EEDGE ) { CPLFree(pszWKT); return false; } if( status != NC_NOERR ) { NCDF_ERR(status); CPLFree(pszWKT); pszWKT = nullptr; } } #ifdef NETCDF_HAS_NC4 else if( m_nWKTNCDFType == NC_STRING ) { char *pszVal = nullptr; int status = nc_get_var1_string(m_nLayerCDFId, m_nWKTVarID, anIndex, &pszVal); if( status == NC_EINVALCOORDS || status == NC_EEDGE ) { return false; } if( status != NC_NOERR ) { NCDF_ERR(status); } else if( pszVal != nullptr ) { pszWKT = CPLStrdup(pszVal); nc_free_string(1, &pszVal); } } #endif if( pszWKT != nullptr ) { OGRGeometry *poGeom = nullptr; CPL_IGNORE_RET_VAL( OGRGeometryFactory::createFromWkt(pszWKT, nullptr, &poGeom)); if( poGeom != nullptr ) { poGeom->assignSpatialReference(GetSpatialRef()); poFeature->SetGeometryDirectly(poGeom); } CPLFree(pszWKT); } } return true; } /************************************************************************/ /* GetNextFeature() */ /************************************************************************/ OGRFeature *netCDFLayer::GetNextFeature() { while( true ) { OGRFeature *poFeature = GetNextRawFeature(); if( poFeature == nullptr ) return nullptr; if( (m_poFilterGeom == nullptr || FilterGeometry(poFeature->GetGeomFieldRef(m_iGeomFieldFilter) )) && (m_poAttrQuery == nullptr || m_poAttrQuery->Evaluate(poFeature)) ) return poFeature; delete poFeature; } } /************************************************************************/ /* GetLayerDefn() */ /************************************************************************/ OGRFeatureDefn *netCDFLayer::GetLayerDefn() { return m_poFeatureDefn; } /************************************************************************/ /* ICreateFeature() */ /************************************************************************/ OGRErr netCDFLayer::ICreateFeature(OGRFeature *poFeature) { m_poDS->SetDefineMode(false); size_t nFeatureIdx = 0; nc_inq_dimlen(m_nLayerCDFId, m_nRecordDimID, &nFeatureIdx); if( m_nProfileDimID >= 0 ) { size_t nProfileCount = 0; nc_inq_dimlen(m_nLayerCDFId, m_nProfileDimID, &nProfileCount); OGRFeature *poProfileToLookup = poFeature->Clone(); poProfileToLookup->SetFID(-1); for( int i = 0; i < m_poFeatureDefn->GetFieldCount(); i++ ) { if( !(poProfileToLookup->IsFieldSetAndNotNull(i)) || m_aoFieldDesc[i].nMainDimId != m_nProfileDimID ) { poProfileToLookup->UnsetField(i); continue; } } OGRGeometry *poGeom = poProfileToLookup->GetGeometryRef(); if( poGeom != nullptr && wkbFlatten(poGeom->getGeometryType()) == wkbPoint ) { poGeom->toPoint()->setZ(0); } size_t nProfileIdx = 0; bool bFoundProfile = false; for( ; nProfileIdx < nProfileCount; nProfileIdx++ ) { int nId = NC_FILL_INT; int status = nc_get_var1_int(m_nLayerCDFId, m_nProfileVarID, &nProfileIdx, &nId); NCDF_ERR(status); if( nId == NC_FILL_INT ) break; OGRFeature *poIterFeature = new OGRFeature(m_poFeatureDefn); if( FillFeatureFromVar(poIterFeature, m_nProfileDimID, nProfileIdx) ) { poGeom = poIterFeature->GetGeometryRef(); if( poGeom != nullptr && wkbFlatten(poGeom->getGeometryType()) == wkbPoint ) { poGeom->toPoint()->setZ(0); } if( poIterFeature->Equal(poProfileToLookup) ) { bFoundProfile = true; delete poIterFeature; break; } } delete poIterFeature; } if( !bFoundProfile ) { if( !m_bProfileVarUnlimited && nProfileIdx == nProfileCount ) { size_t nNewSize = 1 + nProfileCount + nProfileCount / 3; m_poDS->GrowDim(m_nLayerCDFId, m_nProfileDimID, nNewSize); } if( !FillVarFromFeature(poProfileToLookup, m_nProfileDimID, nProfileIdx) ) { delete poProfileToLookup; return OGRERR_FAILURE; } } int nProfileIdIdx = m_poFeatureDefn->GetFieldIndex(m_osProfileDimName); if( nProfileIdIdx < 0 || m_poFeatureDefn->GetFieldDefn(nProfileIdIdx)->GetType() != OFTInteger ) { int nVal = static_cast(nProfileIdx); int status = nc_put_var1_int(m_nLayerCDFId, m_nProfileVarID, &nProfileIdx, &nVal); NCDF_ERR(status); } int nVal = static_cast(nProfileIdx); int status = nc_put_var1_int(m_nLayerCDFId, m_nParentIndexVarID, &nFeatureIdx, &nVal); NCDF_ERR(status); delete poProfileToLookup; } if( !FillVarFromFeature(poFeature, m_nRecordDimID, nFeatureIdx) ) return OGRERR_FAILURE; poFeature->SetFID(nFeatureIdx + 1); return OGRERR_NONE; } /************************************************************************/ /* FillVarFromFeature() */ /************************************************************************/ bool netCDFLayer::FillVarFromFeature(OGRFeature *poFeature, int nMainDimId, size_t nIndex) { size_t anIndex[2]; anIndex[0] = nIndex; anIndex[1] = 0; for( int i = 0; i < m_poFeatureDefn->GetFieldCount(); i++ ) { if( m_aoFieldDesc[i].nMainDimId != nMainDimId ) continue; if( !(poFeature->IsFieldSetAndNotNull(i)) ) { #ifdef NETCDF_HAS_NC4 if( m_bNCDumpCompat && m_aoFieldDesc[i].nType == NC_STRING ) { const char *pszVal = ""; int status = nc_put_var1_string( m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, &pszVal); NCDF_ERR(status); } #endif continue; } int status = NC_NOERR; switch( m_aoFieldDesc[i].nType ) { case NC_CHAR: { const char *pszVal = poFeature->GetFieldAsString(i); if( m_aoFieldDesc[i].nDimCount == 1 ) { if( strlen(pszVal) > 1 && !m_aoFieldDesc[i].bHasWarnedAboutTruncation ) { CPLError(CE_Warning, CPLE_AppDefined, "Content of field %s exceeded the 1 character " "limit and will be truncated", m_poFeatureDefn->GetFieldDefn(i)->GetNameRef()); m_aoFieldDesc[i].bHasWarnedAboutTruncation = true; } status = nc_put_var1_text( m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, pszVal); } else { size_t anCount[2]; anCount[0] = 1; anCount[1] = strlen(pszVal); size_t nWidth = 0; nc_inq_dimlen(m_nLayerCDFId, m_aoFieldDesc[i].nSecDimId, &nWidth); if( anCount[1] > nWidth ) { if( m_bAutoGrowStrings && m_poFeatureDefn->GetFieldDefn(i)->GetWidth() == 0 ) { size_t nNewSize = anCount[1] + anCount[1] / 3; CPLDebug("GDAL_netCDF", "Growing %s from %u to %u", m_poFeatureDefn->GetFieldDefn(i)->GetNameRef(), static_cast(nWidth), static_cast(nNewSize)); m_poDS->GrowDim(m_nLayerCDFId, m_aoFieldDesc[i].nSecDimId, nNewSize); pszVal = poFeature->GetFieldAsString(i); } else { anCount[1] = nWidth; if( !m_aoFieldDesc[i].bHasWarnedAboutTruncation ) { CPLError( CE_Warning, CPLE_AppDefined, "Content of field %s exceeded the %u character " "limit and will be truncated", m_poFeatureDefn->GetFieldDefn(i)->GetNameRef(), static_cast(nWidth)); m_aoFieldDesc[i].bHasWarnedAboutTruncation = true; } } } status = nc_put_vara_text(m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, anCount, pszVal); } break; } #ifdef NETCDF_HAS_NC4 case NC_STRING: { const char *pszVal = poFeature->GetFieldAsString(i); status = nc_put_var1_string(m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, &pszVal); break; } #endif case NC_BYTE: { int nVal = poFeature->GetFieldAsInteger(i); signed char chVal = static_cast(nVal); status = nc_put_var1_schar(m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, &chVal); break; } #ifdef NETCDF_HAS_NC4 case NC_UBYTE: { int nVal = poFeature->GetFieldAsInteger(i); unsigned char uchVal = static_cast(nVal); status = nc_put_var1_uchar(m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, &uchVal); break; } #endif case NC_SHORT: { int nVal = poFeature->GetFieldAsInteger(i); short sVal = static_cast(nVal); status = nc_put_var1_short(m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, &sVal); break; } #ifdef NETCDF_HAS_NC4 case NC_USHORT: { int nVal = poFeature->GetFieldAsInteger(i); unsigned short usVal = static_cast(nVal); status = nc_put_var1_ushort(m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, &usVal); break; } #endif case NC_INT: { int nVal; if( m_poFeatureDefn->GetFieldDefn(i)->GetType() == OFTDate ) { int nYear; int nMonth; int nDay; int nHour; int nMinute; float fSecond; int nTZ; poFeature->GetFieldAsDateTime(i, &nYear, &nMonth, &nDay, &nHour, &nMinute, &fSecond, &nTZ); struct tm brokendowntime; brokendowntime.tm_year = nYear - 1900; brokendowntime.tm_mon = nMonth - 1; brokendowntime.tm_mday = nDay; brokendowntime.tm_hour = 0; brokendowntime.tm_min = 0; brokendowntime.tm_sec = 0; GIntBig nVal64 = CPLYMDHMSToUnixTime(&brokendowntime); if( m_aoFieldDesc[i].bIsDays ) nVal64 /= 86400; nVal = static_cast(nVal64); } else { nVal = poFeature->GetFieldAsInteger(i); } status = nc_put_var1_int(m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, &nVal); break; } #ifdef NETCDF_HAS_NC4 case NC_UINT: { GIntBig nVal = poFeature->GetFieldAsInteger64(i); unsigned int unVal = static_cast(nVal); status = nc_put_var1_uint(m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, &unVal); break; } #endif #ifdef NETCDF_HAS_NC4 case NC_INT64: { GIntBig nVal = poFeature->GetFieldAsInteger64(i); status = nc_put_var1_longlong( m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, &nVal); break; } case NC_UINT64: { double dfVal = poFeature->GetFieldAsDouble(i); GUIntBig nVal = static_cast(dfVal); status = nc_put_var1_ulonglong( m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, &nVal); break; } #endif case NC_FLOAT: { double dfVal = poFeature->GetFieldAsDouble(i); float fVal = static_cast(dfVal); status = nc_put_var1_float(m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, &fVal); break; } case NC_DOUBLE: { double dfVal; if( m_poFeatureDefn->GetFieldDefn(i)->GetType() == OFTDate || m_poFeatureDefn->GetFieldDefn(i)->GetType() == OFTDateTime ) { int nYear; int nMonth; int nDay; int nHour; int nMinute; float fSecond; int nTZ; poFeature->GetFieldAsDateTime(i, &nYear, &nMonth, &nDay, &nHour, &nMinute, &fSecond, &nTZ); struct tm brokendowntime; brokendowntime.tm_year = nYear - 1900; brokendowntime.tm_mon = nMonth - 1; brokendowntime.tm_mday = nDay; brokendowntime.tm_hour = nHour; brokendowntime.tm_min = nMinute; brokendowntime.tm_sec = static_cast(fSecond); GIntBig nVal = CPLYMDHMSToUnixTime(&brokendowntime); dfVal = static_cast(nVal) + fmod(fSecond, 1.0f); if( m_aoFieldDesc[i].bIsDays ) dfVal /= 86400.0; } else { dfVal = poFeature->GetFieldAsDouble(i); } status = nc_put_var1_double(m_nLayerCDFId, m_aoFieldDesc[i].nVarId, anIndex, &dfVal); break; } default: break; } NCDF_ERR(status); if( status != NC_NOERR ) { return false; } } OGRGeometry *poGeom = poFeature->GetGeometryRef(); if( wkbFlatten(m_poFeatureDefn->GetGeomType()) == wkbPoint && poGeom != nullptr && wkbFlatten(poGeom->getGeometryType()) == wkbPoint ) { if( m_osProfileDimName.empty() || nMainDimId == m_nProfileDimID ) { auto poPoint = poGeom->toPoint(); double dfX = poPoint->getX(); double dfY = poPoint->getY(); int status; if( m_nXVarNCDFType == NC_DOUBLE ) status = nc_put_var1_double(m_nLayerCDFId, m_nXVarID, anIndex, &dfX); else { float fX = static_cast(dfX); status = nc_put_var1_float(m_nLayerCDFId, m_nXVarID, anIndex, &fX); } NCDF_ERR(status); if( status != NC_NOERR ) { return false; } if( m_nYVarNCDFType == NC_DOUBLE ) status = nc_put_var1_double(m_nLayerCDFId, m_nYVarID, anIndex, &dfY); else { float fY = static_cast(dfY); status = nc_put_var1_float(m_nLayerCDFId, m_nYVarID, anIndex, &fY); } NCDF_ERR(status); if( status != NC_NOERR ) { return false; } } if( m_poFeatureDefn->GetGeomType() == wkbPoint25D && (m_osProfileDimName.empty() || nMainDimId == m_nRecordDimID) ) { int status; double dfZ = poGeom->toPoint()->getZ(); if( m_nZVarNCDFType == NC_DOUBLE ) status = nc_put_var1_double(m_nLayerCDFId, m_nZVarID, anIndex, &dfZ); else { float fZ = static_cast(dfZ); status = nc_put_var1_float(m_nLayerCDFId, m_nZVarID, anIndex, &fZ); } NCDF_ERR(status); if( status != NC_NOERR ) { return false; } } } else if( m_poFeatureDefn->GetGeomType() != wkbNone && m_nWKTVarID >= 0 && poGeom != nullptr ) { char *pszWKT = nullptr; poGeom->exportToWkt(&pszWKT, wkbVariantIso); int status; #ifdef NETCDF_HAS_NC4 if( m_nWKTNCDFType == NC_STRING ) { const char *pszWKTConst = pszWKT; status = nc_put_var1_string(m_nLayerCDFId, m_nWKTVarID, anIndex, &pszWKTConst); } else #endif { size_t anCount[2]; anCount[0] = 1; anCount[1] = strlen(pszWKT); if( anCount[1] > static_cast(m_nWKTMaxWidth) ) { if( m_bAutoGrowStrings ) { size_t nNewSize = anCount[1] + anCount[1] / 3; CPLDebug("GDAL_netCDF", "Growing %s from %u to %u", m_osWKTVarName.c_str(), static_cast(m_nWKTMaxWidth), static_cast(nNewSize)); m_poDS->GrowDim(m_nLayerCDFId, m_nWKTMaxWidthDimId, nNewSize); m_nWKTMaxWidth = static_cast(nNewSize); status = nc_put_vara_text(m_nLayerCDFId, m_nWKTVarID, anIndex, anCount, pszWKT); } else { CPLError(CE_Failure, CPLE_AppDefined, "Cannot write geometry as WKT. Would require %d " "characters but field width is %d", static_cast(anCount[1]), m_nWKTMaxWidth); status = NC_NOERR; } } else { status = nc_put_vara_text(m_nLayerCDFId, m_nWKTVarID, anIndex, anCount, pszWKT); } } CPLFree(pszWKT); NCDF_ERR(status); if( status != NC_NOERR ) { return false; } } #ifdef NETCDF_HAS_NC4 else if( m_poFeatureDefn->GetGeomType() != wkbNone && m_nWKTVarID >= 0 && poGeom == nullptr && m_nWKTNCDFType == NC_STRING && m_bNCDumpCompat ) { const char *pszWKTConst = ""; int status = nc_put_var1_string(m_nLayerCDFId, m_nWKTVarID, anIndex, &pszWKTConst); NCDF_ERR(status); } #endif return true; } /************************************************************************/ /* AddField() */ /************************************************************************/ bool netCDFLayer::AddField(int nVarID) { if( nVarID == m_nWKTVarID ) return false; char szName[NC_MAX_NAME + 1]; szName[0] = '\0'; CPL_IGNORE_RET_VAL(nc_inq_varname(m_nLayerCDFId, nVarID, szName)); nc_type vartype = NC_NAT; nc_inq_vartype(m_nLayerCDFId, nVarID, &vartype); OGRFieldType eType = OFTString; OGRFieldSubType eSubType = OFSTNone; int nWidth = 0; NCDFNoDataUnion nodata; memset(&nodata, 0, sizeof(nodata)); int nDimCount = 1; nc_inq_varndims(m_nLayerCDFId, nVarID, &nDimCount); int anDimIds[2] = { -1, -1 }; if( (vartype == NC_CHAR && nDimCount <= 2) || (vartype != NC_CHAR && nDimCount == 1) ) { nc_inq_vardimid(m_nLayerCDFId, nVarID, anDimIds); } else { return false; } switch( vartype ) { case NC_BYTE: { eType = OFTInteger; char *pszValue = nullptr; if( GetFillValue(nVarID, &pszValue) == CE_None ) nodata.chVal = static_cast(atoi(pszValue)); else nodata.chVal = NC_FILL_BYTE; CPLFree(pszValue); break; } #ifdef NETCDF_HAS_NC4 case NC_UBYTE: { eType = OFTInteger; char *pszValue = nullptr; if( GetFillValue(nVarID, &pszValue) == CE_None ) nodata.uchVal = static_cast(atoi(pszValue)); else nodata.uchVal = NC_FILL_UBYTE; CPLFree(pszValue); break; } #endif case NC_CHAR: { eType = OFTString; if( nDimCount == 1 ) { nWidth = 1; } else if( nDimCount == 2 ) { size_t nDimLen = 0; nc_inq_dimlen(m_nLayerCDFId, anDimIds[1], &nDimLen); nWidth = static_cast(nDimLen); } break; } #ifdef NETCDF_HAS_NC4 case NC_STRING: { eType = OFTString; break; } #endif case NC_SHORT: { eType = OFTInteger; eSubType = OFSTInt16; char *pszValue = nullptr; if( GetFillValue(nVarID, &pszValue) == CE_None ) nodata.sVal = static_cast(atoi(pszValue)); else nodata.sVal = NC_FILL_SHORT; CPLFree(pszValue); break; } #ifdef NETCDF_HAS_NC4 case NC_USHORT: { eType = OFTInteger; char *pszValue = nullptr; if( GetFillValue(nVarID, &pszValue) == CE_None ) nodata.usVal = static_cast(atoi(pszValue)); else nodata.usVal = NC_FILL_USHORT; CPLFree(pszValue); break; } #endif case NC_INT: { eType = OFTInteger; char *pszValue = nullptr; if( GetFillValue(nVarID, &pszValue) == CE_None ) nodata.nVal = atoi(pszValue); else nodata.nVal = NC_FILL_INT; CPLFree(pszValue); break; } #ifdef NETCDF_HAS_NC4 case NC_UINT: { eType = OFTInteger64; char *pszValue = nullptr; if( GetFillValue(nVarID, &pszValue) == CE_None ) nodata.unVal = static_cast(CPLAtoGIntBig(pszValue)); else nodata.unVal = NC_FILL_UINT; CPLFree(pszValue); break; } #endif #ifdef NETCDF_HAS_NC4 case NC_INT64: { eType = OFTInteger64; char *pszValue = nullptr; if( GetFillValue(nVarID, &pszValue) == CE_None ) nodata.nVal64 = CPLAtoGIntBig(pszValue); else nodata.nVal64 = NC_FILL_INT64; CPLFree(pszValue); break; } case NC_UINT64: { eType = OFTReal; char *pszValue = nullptr; if( GetFillValue(nVarID, &pszValue) == CE_None ) { nodata.unVal64 = 0; for( int i = 0; pszValue[i] != '\0'; i++ ) { nodata.unVal64 = nodata.unVal64 * 10 + (pszValue[i] - '0'); } } else nodata.unVal64 = NC_FILL_UINT64; CPLFree(pszValue); break; } #endif case NC_FLOAT: { eType = OFTReal; eSubType = OFSTFloat32; double dfValue; if( GetFillValue(nVarID, &dfValue) == CE_None ) nodata.fVal = static_cast(dfValue); else nodata.fVal = NC_FILL_FLOAT; break; } case NC_DOUBLE: { eType = OFTReal; double dfValue; if( GetFillValue(nVarID, &dfValue) == CE_None ) nodata.dfVal = dfValue; else nodata.dfVal = NC_FILL_DOUBLE; break; } default: { CPLDebug("GDAL_netCDF", "Variable %s has type %d, which is unhandled", szName, vartype); return false; } } bool bIsDays = false; char *pszValue = nullptr; if( NCDFGetAttr(m_nLayerCDFId, nVarID, "ogr_field_type", &pszValue) == CE_None ) { if( (eType == OFTInteger || eType == OFTReal) && EQUAL(pszValue, "Date") ) { eType = OFTDate; // cppcheck-suppress knownConditionTrueFalse bIsDays = (eType == OFTInteger); } else if( (eType == OFTInteger || eType == OFTReal) && EQUAL(pszValue, "DateTime") ) eType = OFTDateTime; else if( eType == OFTReal && EQUAL(pszValue, "Integer64") ) eType = OFTInteger64; else if( eType == OFTInteger && EQUAL(pszValue, "Integer(Boolean)") ) eSubType = OFSTBoolean; } CPLFree(pszValue); pszValue = nullptr; if( NCDFGetAttr(m_nLayerCDFId, nVarID, "units", &pszValue) == CE_None ) { if( (eType == OFTInteger || eType == OFTReal || eType == OFTDate) && (EQUAL(pszValue, "seconds since 1970-1-1 0:0:0") || EQUAL(pszValue, "seconds since 1970-01-01 00:00:00")) ) { if( eType != OFTDate ) eType = OFTDateTime; bIsDays = false; } else if( (eType == OFTInteger || eType == OFTReal || eType == OFTDate) && (EQUAL(pszValue, "days since 1970-1-1") || EQUAL(pszValue, "days since 1970-01-01")) ) { eType = OFTDate; bIsDays = true; } } CPLFree(pszValue); pszValue = nullptr; if( NCDFGetAttr(m_nLayerCDFId, nVarID, "ogr_field_name", &pszValue) == CE_None ) { snprintf(szName, sizeof(szName), "%s", pszValue); } CPLFree(pszValue); pszValue = nullptr; if( NCDFGetAttr(m_nLayerCDFId, nVarID, "ogr_field_width", &pszValue) == CE_None ) { nWidth = atoi(pszValue); } CPLFree(pszValue); pszValue = nullptr; int nPrecision = 0; if( NCDFGetAttr(m_nLayerCDFId, nVarID, "ogr_field_precision", &pszValue) == CE_None ) { nPrecision = atoi(pszValue); } CPLFree(pszValue); /* pszValue = NULL; */ OGRFieldDefn oFieldDefn(szName, eType); oFieldDefn.SetSubType(eSubType); oFieldDefn.SetWidth(nWidth); oFieldDefn.SetPrecision(nPrecision); FieldDesc fieldDesc; fieldDesc.uNoData = nodata; fieldDesc.nType = vartype; fieldDesc.nVarId = nVarID; fieldDesc.nDimCount = nDimCount; fieldDesc.nMainDimId = anDimIds[0]; fieldDesc.nSecDimId = anDimIds[1]; fieldDesc.bHasWarnedAboutTruncation = false; fieldDesc.bIsDays = bIsDays; m_aoFieldDesc.push_back(fieldDesc); m_poFeatureDefn->AddFieldDefn(&oFieldDefn); return true; } /************************************************************************/ /* CreateField() */ /************************************************************************/ OGRErr netCDFLayer::CreateField(OGRFieldDefn *poFieldDefn, int /* bApproxOK */) { int nSecDimId = -1; int nVarID = -1; int status; const netCDFWriterConfigField *poConfig = nullptr; if( m_poDS->oWriterConfig.m_bIsValid ) { std::map::const_iterator oIter; if( m_poLayerConfig != nullptr && (oIter = m_poLayerConfig->m_oFields.find( poFieldDefn->GetNameRef())) != m_poLayerConfig->m_oFields.end() ) { poConfig = &(oIter->second); } else if( (oIter = m_poDS->oWriterConfig.m_oFields.find( poFieldDefn->GetNameRef())) != m_poDS->oWriterConfig.m_oFields.end() ) { poConfig = &(oIter->second); } } if( !m_osProfileDimName.empty() && EQUAL(poFieldDefn->GetNameRef(), m_osProfileDimName) && poFieldDefn->GetType() == OFTInteger ) { FieldDesc fieldDesc; fieldDesc.uNoData.nVal = NC_FILL_INT; fieldDesc.nType = NC_INT; fieldDesc.nVarId = m_nProfileVarID; fieldDesc.nDimCount = 1; fieldDesc.nMainDimId = m_nProfileDimID; fieldDesc.nSecDimId = -1; fieldDesc.bHasWarnedAboutTruncation = false; fieldDesc.bIsDays = false; m_aoFieldDesc.push_back(fieldDesc); m_poFeatureDefn->AddFieldDefn(poFieldDefn); return OGRERR_NONE; } m_poDS->SetDefineMode(true); // Try to use the field name as variable name, but detects conflict first CPLString osVarName(poConfig != nullptr ? poConfig->m_osNetCDFName : CPLString(poFieldDefn->GetNameRef())); status = nc_inq_varid(m_nLayerCDFId, osVarName, &nVarID); if( status == NC_NOERR ) { for( int i = 1; i <= 100; i++ ) { osVarName = CPLSPrintf("%s%d", poFieldDefn->GetNameRef(), i); status = nc_inq_varid(m_nLayerCDFId, osVarName, &nVarID); if( status != NC_NOERR ) break; } CPLDebug("netCDF", "Field %s is written in variable %s", poFieldDefn->GetNameRef(), osVarName.c_str()); } const char *pszVarName = osVarName.c_str(); NCDFNoDataUnion nodata; memset(&nodata, 0, sizeof(nodata)); const OGRFieldType eType = poFieldDefn->GetType(); const OGRFieldSubType eSubType = poFieldDefn->GetSubType(); nc_type nType = NC_NAT; int nDimCount = 1; // Find which is the dimension that this variable should be indexed against int nMainDimId = m_nRecordDimID; if( !m_osProfileVariables.empty() ) { char **papszTokens = CSLTokenizeString2(m_osProfileVariables, ",", CSLT_HONOURSTRINGS); if( CSLFindString(papszTokens, poFieldDefn->GetNameRef()) >= 0 ) nMainDimId = m_nProfileDimID; CSLDestroy(papszTokens); } if( poConfig != nullptr && !poConfig->m_osMainDim.empty() ) { int ndims = 0; status = nc_inq_ndims(m_nLayerCDFId, &ndims); NCDF_ERR(status); bool bFound = false; for( int idim = 0; idim < ndims; idim++ ) { char szDimName[NC_MAX_NAME + 1]; szDimName[0] = 0; status = nc_inq_dimname(m_poDS->cdfid, idim, szDimName); NCDF_ERR(status); if( strcmp(poConfig->m_osMainDim, szDimName) == 0 ) { nMainDimId = idim; bFound = true; break; } } if( !bFound ) { CPLError(CE_Failure, CPLE_AppDefined, "Dimension '%s' does not exist", poConfig->m_osMainDim.c_str()); } } switch( eType ) { case OFTString: case OFTStringList: case OFTIntegerList: case OFTRealList: { if( poFieldDefn->GetWidth() == 1 ) { nType = NC_CHAR; status = nc_def_var(m_nLayerCDFId, pszVarName, nType, 1, &nMainDimId, &nVarID); } #ifdef NETCDF_HAS_NC4 else if( m_poDS->eFormat == NCDF_FORMAT_NC4 && m_bUseStringInNC4 ) { nType = NC_STRING; status = nc_def_var(m_nLayerCDFId, pszVarName, nType, 1, &nMainDimId, &nVarID); } #endif else { if( poFieldDefn->GetWidth() == 0 && !m_bAutoGrowStrings ) { if( m_nDefaultMaxWidthDimId < 0 ) { status = nc_def_dim(m_nLayerCDFId, "string_default_max_width", m_nDefaultWidth, &m_nDefaultMaxWidthDimId); NCDF_ERR(status); if( status != NC_NOERR ) { return OGRERR_FAILURE; } } nSecDimId = m_nDefaultMaxWidthDimId; } else { size_t nDim = poFieldDefn->GetWidth() == 0 ? m_nDefaultWidth : poFieldDefn->GetWidth(); status = nc_def_dim(m_nLayerCDFId, CPLSPrintf("%s_max_width", pszVarName), nDim, &nSecDimId); NCDF_ERR(status); if( status != NC_NOERR ) { return OGRERR_FAILURE; } } nDimCount = 2; int anDims[2] = {nMainDimId, nSecDimId}; nType = NC_CHAR; status = nc_def_var(m_nLayerCDFId, pszVarName, nType, 2, anDims, &nVarID); } NCDF_ERR(status); if( status != NC_NOERR ) { return OGRERR_FAILURE; } break; } case OFTInteger: { nType = eSubType == OFSTBoolean ? NC_BYTE : (eSubType == OFSTInt16) ? NC_SHORT : NC_INT; if( nType == NC_BYTE ) nodata.chVal = NC_FILL_BYTE; else if( nType == NC_SHORT ) nodata.sVal = NC_FILL_SHORT; else if( nType == NC_INT ) nodata.nVal = NC_FILL_INT; status = nc_def_var(m_nLayerCDFId, pszVarName, nType, 1, &nMainDimId, &nVarID); NCDF_ERR(status); if( status != NC_NOERR ) { return OGRERR_FAILURE; } if( eSubType == OFSTBoolean ) { signed char anRange[2] = { 0, 1 }; nc_put_att_schar(m_nLayerCDFId, nVarID, "valid_range", NC_BYTE, 2, anRange); } break; } case OFTInteger64: { nType = NC_DOUBLE; nodata.dfVal = NC_FILL_DOUBLE; #ifdef NETCDF_HAS_NC4 if( m_poDS->eFormat == NCDF_FORMAT_NC4 ) { nType = NC_INT64; nodata.nVal64 = NC_FILL_INT64; } #endif status = nc_def_var(m_nLayerCDFId, pszVarName, nType, 1, &nMainDimId, &nVarID); NCDF_ERR(status); if( status != NC_NOERR ) { return OGRERR_FAILURE; } break; } case OFTReal: { nType = (eSubType == OFSTFloat32) ? NC_FLOAT : NC_DOUBLE; if( eSubType == OFSTFloat32 ) nodata.fVal = NC_FILL_FLOAT; else nodata.dfVal = NC_FILL_DOUBLE; status = nc_def_var(m_nLayerCDFId, pszVarName, nType, 1, &nMainDimId, &nVarID); NCDF_ERR(status); if( status != NC_NOERR ) { return OGRERR_FAILURE; } break; } case OFTDate: { nType = NC_INT; status = nc_def_var(m_nLayerCDFId, pszVarName, nType, 1, &nMainDimId, &nVarID); NCDF_ERR(status); if( status != NC_NOERR ) { return OGRERR_FAILURE; } nodata.nVal = NC_FILL_INT; status = nc_put_att_text(m_nLayerCDFId, nVarID, CF_UNITS, strlen("days since 1970-1-1"), "days since 1970-1-1"); NCDF_ERR(status); break; } case OFTDateTime: { nType = NC_DOUBLE; status = nc_def_var(m_nLayerCDFId, pszVarName, nType, 1, &nMainDimId, &nVarID); NCDF_ERR(status); if( status != NC_NOERR ) { return OGRERR_FAILURE; } nodata.dfVal = NC_FILL_DOUBLE; status = nc_put_att_text(m_nLayerCDFId, nVarID, CF_UNITS, strlen("seconds since 1970-1-1 0:0:0"), "seconds since 1970-1-1 0:0:0"); NCDF_ERR(status); break; } default: return OGRERR_FAILURE; } FieldDesc fieldDesc; fieldDesc.uNoData = nodata; fieldDesc.nType = nType; fieldDesc.nVarId = nVarID; fieldDesc.nDimCount = nDimCount; fieldDesc.nMainDimId = nMainDimId; fieldDesc.nSecDimId = nSecDimId; fieldDesc.bHasWarnedAboutTruncation = false; fieldDesc.bIsDays = (eType == OFTDate); m_aoFieldDesc.push_back(fieldDesc); const char *pszLongName = CPLSPrintf("Field %s", poFieldDefn->GetNameRef()); status = nc_put_att_text(m_nLayerCDFId, nVarID, CF_LNG_NAME, strlen(pszLongName), pszLongName); NCDF_ERR(status); if( m_bWriteGDALTags ) { status = nc_put_att_text(m_nLayerCDFId, nVarID, "ogr_field_name", strlen(poFieldDefn->GetNameRef()), poFieldDefn->GetNameRef()); NCDF_ERR(status); const char *pszType = OGRFieldDefn::GetFieldTypeName(eType); if( eSubType != OFSTNone ) { pszType = CPLSPrintf("%s(%s)", pszType, OGRFieldDefn::GetFieldSubTypeName(eSubType)); } status = nc_put_att_text(m_nLayerCDFId, nVarID, "ogr_field_type", strlen(pszType), pszType); NCDF_ERR(status); const int nWidth = poFieldDefn->GetWidth(); if(nWidth || nType == NC_CHAR) { status = nc_put_att_int(m_nLayerCDFId, nVarID, "ogr_field_width", NC_INT, 1, &nWidth); NCDF_ERR(status); const int nPrecision = poFieldDefn->GetPrecision(); if( nPrecision ) { status = nc_put_att_int(m_nLayerCDFId, nVarID, "ogr_field_precision", NC_INT, 1, &nPrecision); NCDF_ERR(status); } } } // nc_put_att_text(m_nLayerCDFId, nVarID, CF_UNITS, // strlen("none"), "none"); if( !m_osGridMapping.empty() && nMainDimId == m_nRecordDimID ) { status = nc_put_att_text(m_nLayerCDFId, nVarID, CF_GRD_MAPPING, m_osGridMapping.size(), m_osGridMapping.c_str()); NCDF_ERR(status); } if( !m_osCoordinatesValue.empty() && nMainDimId == m_nRecordDimID ) { status = nc_put_att_text(m_nLayerCDFId, nVarID, CF_COORDINATES, m_osCoordinatesValue.size(), m_osCoordinatesValue.c_str()); NCDF_ERR(status); } if( poConfig != nullptr ) { netCDFWriteAttributesFromConf(m_nLayerCDFId, nVarID, poConfig->m_aoAttributes); } m_poFeatureDefn->AddFieldDefn(poFieldDefn); return OGRERR_NONE; } /************************************************************************/ /* GetFeatureCount() */ /************************************************************************/ GIntBig netCDFLayer::GetFeatureCount(int bForce) { if( m_poFilterGeom == nullptr && m_poAttrQuery == nullptr ) { size_t nDimLen; nc_inq_dimlen(m_nLayerCDFId, m_nRecordDimID, &nDimLen); return static_cast(nDimLen); } return OGRLayer::GetFeatureCount(bForce); } /************************************************************************/ /* TestCapability() */ /************************************************************************/ int netCDFLayer::TestCapability(const char *pszCap) { if( EQUAL(pszCap, OLCSequentialWrite) ) return m_poDS->GetAccess() == GA_Update; if( EQUAL(pszCap, OLCCreateField) ) return m_poDS->GetAccess() == GA_Update; if( EQUAL(pszCap, OLCFastFeatureCount) ) return m_poFilterGeom == nullptr && m_poAttrQuery == nullptr; return FALSE; }