1 /******************************************************************************
2  *
3  * Project:  netCDF read/write Driver
4  * Author:   Even Rouault <even.rouault at spatialys.com>
5  *
6  ******************************************************************************
7  * Copyright (c) 2019, Even Rouault <even.rouault at spatialys.com>
8  *
9  * Permission is hereby granted, free of charge, to any person obtaining a
10  * copy of this software and associated documentation files (the "Software"),
11  * to deal in the Software without restriction, including without limitation
12  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
13  * and/or sell copies of the Software, and to permit persons to whom the
14  * Software is furnished to do so, subject to the following conditions:
15  *
16  * The above copyright notice and this permission notice shall be included
17  * in all copies or substantial portions of the Software.
18  *
19  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
20  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
22  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
24  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
25  * DEALINGS IN THE SOFTWARE.
26  ****************************************************************************/
27 
28 #include <algorithm>
29 #include <limits>
30 #include <map>
31 
32 #include "netcdfdataset.h"
33 
34 #ifdef NETCDF_HAS_NC4
35 
36 /************************************************************************/
37 /*                         netCDFSharedResources                        */
38 /************************************************************************/
39 
40 class netCDFSharedResources
41 {
42     friend class netCDFDataset;
43 
44     bool m_bImappIsInElements = true;
45     bool m_bReadOnly = true;
46     int m_cdfid = 0;
47 #ifdef ENABLE_NCDUMP
48     bool m_bFileToDestroyAtClosing = false;
49 #endif
50     CPLString m_osFilename{};
51 #ifdef ENABLE_UFFD
52     cpl_uffd_context *m_pUffdCtx = nullptr;
53 #endif
54     bool m_bDefineMode = false;
55     std::map<int, int> m_oMapDimIdToGroupId{};
56     bool m_bIsInIndexingVariable = false;
57 
58 public:
59     netCDFSharedResources();
60     ~netCDFSharedResources();
61 
GetCDFId() const62     inline int GetCDFId() const { return m_cdfid; }
IsReadOnly() const63     inline bool IsReadOnly() const { return m_bReadOnly; }
64     bool SetDefineMode(bool bNewDefineMode);
65     int GetBelongingGroupOfDim(int startgid, int dimid);
GetImappIsInElements() const66     inline bool GetImappIsInElements() const { return m_bImappIsInElements; }
67 
SetIsInGetIndexingVariable(bool b)68     void SetIsInGetIndexingVariable(bool b) { m_bIsInIndexingVariable = b; }
GetIsInIndexingVariable() const69     bool GetIsInIndexingVariable() const { return m_bIsInIndexingVariable; }
70 };
71 
72 /************************************************************************/
73 /*                       netCDFSharedResources()                        */
74 /************************************************************************/
75 
netCDFSharedResources()76 netCDFSharedResources::netCDFSharedResources()
77 {
78     // netcdf >= 4.4 uses imapp argument of nc_get/put_varm as a stride in
79     // elements, whereas earlier versions use bytes.
80     CPLStringList aosVersionNumbers(CSLTokenizeString2(nc_inq_libvers(), ".", 0));
81     m_bImappIsInElements = false;
82     if( aosVersionNumbers.size() >= 3 )
83     {
84         m_bImappIsInElements =
85             (atoi(aosVersionNumbers[0]) > 4 || atoi(aosVersionNumbers[1]) >= 4);
86     }
87 }
88 
89 /************************************************************************/
90 /*                       GetBelongingGroupOfDim()                       */
91 /************************************************************************/
92 
GetBelongingGroupOfDim(int startgid,int dimid)93 int netCDFSharedResources::GetBelongingGroupOfDim(int startgid, int dimid)
94 {
95     // Am I missing a netCDF API to do this directly ?
96     auto oIter = m_oMapDimIdToGroupId.find(dimid);
97     if (oIter != m_oMapDimIdToGroupId.end() )
98         return oIter->second;
99 
100     int gid = startgid;
101     while( true )
102     {
103         int nbDims = 0;
104         NCDF_ERR(nc_inq_ndims(gid, &nbDims));
105         if( nbDims > 0 )
106         {
107             std::vector<int> dimids(nbDims);
108             NCDF_ERR(nc_inq_dimids(gid, &nbDims, &dimids[0], FALSE));
109             for( int i = 0; i < nbDims; i++ )
110             {
111                 m_oMapDimIdToGroupId[dimid] = gid;
112                 if( dimids[i] == dimid )
113                     return gid;
114             }
115         }
116         int nParentGID = 0;
117         if( nc_inq_grp_parent(gid, &nParentGID) != NC_NOERR )
118             return startgid;
119         gid = nParentGID;
120     }
121 }
122 
123 
124 /************************************************************************/
125 /*                            SetDefineMode()                           */
126 /************************************************************************/
127 
SetDefineMode(bool bNewDefineMode)128 bool netCDFSharedResources::SetDefineMode( bool bNewDefineMode )
129 {
130     // Do nothing if already in new define mode
131     // or if dataset is in read-only mode.
132     if( m_bDefineMode == bNewDefineMode || m_bReadOnly )
133         return true;
134 
135     CPLDebug("GDAL_netCDF", "SetDefineMode(%d) old=%d",
136              static_cast<int>(bNewDefineMode), static_cast<int>(m_bDefineMode));
137 
138     m_bDefineMode = bNewDefineMode;
139 
140     int status;
141     if( m_bDefineMode )
142         status = nc_redef(m_cdfid);
143     else
144         status = nc_enddef(m_cdfid);
145 
146     NCDF_ERR(status);
147     return status == NC_NOERR;
148 }
149 
150 /************************************************************************/
151 /*                           netCDFGroup                                */
152 /************************************************************************/
153 
154 class netCDFGroup final: public GDALGroup
155 {
156     std::shared_ptr<netCDFSharedResources> m_poShared;
157     int m_gid = 0;
158     CPLStringList m_aosStructuralInfo{};
159 
retrieveName(int gid)160     static std::string retrieveName(int gid)
161     {
162         CPLMutexHolderD(&hNCMutex);
163         char szName[NC_MAX_NAME + 1] = {};
164         NCDF_ERR(nc_inq_grpname(gid, szName));
165         return szName;
166     }
167 
168 public:
169     netCDFGroup(const std::shared_ptr<netCDFSharedResources>& poShared, int gid);
170 
171     std::vector<std::string> GetGroupNames(CSLConstList papszOptions) const override;
172     std::shared_ptr<GDALGroup> OpenGroup(const std::string& osName,
173                                          CSLConstList papszOptions) const override;
174 
175     std::vector<std::string> GetMDArrayNames(CSLConstList papszOptions) const override;
176     std::shared_ptr<GDALMDArray> OpenMDArray(const std::string& osName,
177                                              CSLConstList papszOptions) const override;
178 
179     std::vector<std::shared_ptr<GDALDimension>> GetDimensions(CSLConstList papszOptions) const override;
180 
181     std::shared_ptr<GDALAttribute> GetAttribute(const std::string& osName) const override;
182 
183     std::vector<std::shared_ptr<GDALAttribute>> GetAttributes(CSLConstList papszOptions) const override;
184 
185     std::shared_ptr<GDALGroup> CreateGroup(const std::string& osName,
186                                             CSLConstList papszOptions) override;
187 
188     std::shared_ptr<GDALDimension> CreateDimension(const std::string& osName,
189                                                     const std::string& osType,
190                                                     const std::string& osDirection,
191                                                     GUInt64 nSize,
192                                                     CSLConstList papszOptions) override;
193 
194     std::shared_ptr<GDALMDArray> CreateMDArray(const std::string& osName,
195                                                const std::vector<std::shared_ptr<GDALDimension>>& aoDimensions,
196                                                const GDALExtendedDataType& oDataType,
197                                                CSLConstList papszOptions) override;
198 
199     std::shared_ptr<GDALAttribute> CreateAttribute(
200         const std::string& osName,
201         const std::vector<GUInt64>& anDimensions,
202         const GDALExtendedDataType& oDataType,
203         CSLConstList papszOptions) override;
204 
205     CSLConstList GetStructuralInfo() const override;
206 };
207 
208 /************************************************************************/
209 /*                   netCDFVirtualGroupBySameDimension                  */
210 /************************************************************************/
211 
212 class netCDFVirtualGroupBySameDimension final: public GDALGroup
213 {
214     // the real group to which we derived this virtual group from
215     std::shared_ptr<netCDFGroup> m_poGroup;
216     std::string m_osDimName{};
217 
218 public:
219     netCDFVirtualGroupBySameDimension(const std::shared_ptr<netCDFGroup>& poGroup,
220                                       const std::string& osDimName);
221 
222     std::vector<std::string> GetMDArrayNames(CSLConstList papszOptions) const override;
223     std::shared_ptr<GDALMDArray> OpenMDArray(const std::string& osName,
224                                              CSLConstList papszOptions) const override;
225 };
226 
227 /************************************************************************/
228 /*                         netCDFDimension                              */
229 /************************************************************************/
230 
231 class netCDFDimension final: public GDALDimension
232 {
233     std::shared_ptr<netCDFSharedResources> m_poShared;
234     int m_gid = 0;
235     int m_dimid = 0;
236 
retrieveName(int cfid,int dimid)237     static std::string retrieveName(int cfid, int dimid)
238     {
239         CPLMutexHolderD(&hNCMutex);
240         char szName[NC_MAX_NAME + 1] = {};
241         NCDF_ERR(nc_inq_dimname(cfid, dimid, szName));
242         return szName;
243     }
244 
retrieveSize(int cfid,int dimid)245     static GUInt64 retrieveSize(int cfid, int dimid)
246     {
247         CPLMutexHolderD(&hNCMutex);
248         size_t nDimLen = 0;
249         NCDF_ERR(nc_inq_dimlen(cfid, dimid, &nDimLen));
250         return nDimLen;
251     }
252 
253 public:
254     netCDFDimension(const std::shared_ptr<netCDFSharedResources>& poShared, int cfid,
255                     int dimid, size_t nForcedSize, const std::string& osType);
256 
257     std::shared_ptr<GDALMDArray> GetIndexingVariable() const override;
258 
GetId() const259     int GetId() const { return m_dimid; }
260 };
261 
262 /************************************************************************/
263 /*                         netCDFAttribute                              */
264 /************************************************************************/
265 
266 
267 class netCDFAttribute final: public GDALAttribute
268 {
269     std::shared_ptr<netCDFSharedResources> m_poShared;
270     int m_gid = 0;
271     int m_varid = 0;
272     size_t m_nTextLength = 0;
273     std::vector<std::shared_ptr<GDALDimension>> m_dims{};
274     nc_type m_nAttType = NC_NAT;
275     mutable std::unique_ptr<GDALExtendedDataType> m_dt;
276     mutable bool m_bPerfectDataTypeMatch = false;
277 
278 protected:
279     netCDFAttribute(const std::shared_ptr<netCDFSharedResources>& poShared,
280                     int gid, int varid, const std::string& name);
281 
282     netCDFAttribute(const std::shared_ptr<netCDFSharedResources>& poShared,
283                     int gid, int varid, const std::string& osName,
284                     const std::vector<GUInt64>& anDimensions,
285                     const GDALExtendedDataType& oDataType,
286                     CSLConstList papszOptions);
287 
288     bool IRead(const GUInt64* arrayStartIdx,     // array of size GetDimensionCount()
289                       const size_t* count,                 // array of size GetDimensionCount()
290                       const GInt64* arrayStep,        // step in elements
291                       const GPtrDiff_t* bufferStride, // stride in elements
292                       const GDALExtendedDataType& bufferDataType,
293                       void* pDstBuffer) const override;
294 
295     bool IWrite(const GUInt64* arrayStartIdx,     // array of size GetDimensionCount()
296                       const size_t* count,                 // array of size GetDimensionCount()
297                       const GInt64* arrayStep,        // step in elements
298                       const GPtrDiff_t* bufferStride, // stride in elements
299                       const GDALExtendedDataType& bufferDataType,
300                       const void* pSrcBuffer) override;
301 
302 public:
303     static std::shared_ptr<netCDFAttribute> Create(
304                     const std::shared_ptr<netCDFSharedResources>& poShared,
305                     int gid, int varid, const std::string& name);
306 
307     static std::shared_ptr<netCDFAttribute> Create(
308                     const std::shared_ptr<netCDFSharedResources>& poShared,
309                     int gid, int varid, const std::string& osName,
310                     const std::vector<GUInt64>& anDimensions,
311                     const GDALExtendedDataType& oDataType,
312                     CSLConstList papszOptions);
313 
GetDimensions() const314     const std::vector<std::shared_ptr<GDALDimension>>& GetDimensions() const override { return m_dims; }
315 
316     const GDALExtendedDataType &GetDataType() const override;
317 };
318 
319 /************************************************************************/
320 /*                         netCDFVariable                               */
321 /************************************************************************/
322 
323 class netCDFVariable final: public GDALMDArray
324 {
325     std::shared_ptr<netCDFSharedResources> m_poShared;
326     int m_gid = 0;
327     int m_varid = 0;
328     int m_nDims = 0;
329     mutable std::vector<std::shared_ptr<GDALDimension>> m_dims{};
330     mutable nc_type m_nVarType = NC_NAT;
331     mutable std::unique_ptr<GDALExtendedDataType> m_dt;
332     mutable bool m_bPerfectDataTypeMatch = false;
333     mutable std::vector<GByte> m_abyNoData{};
334     mutable bool m_bGetRawNoDataValueHasRun = false;
335     bool m_bHasWrittenData = true;
336     std::string m_osUnit{};
337     CPLStringList m_aosStructuralInfo{};
338     mutable bool m_bSRSRead = false;
339     mutable std::shared_ptr<OGRSpatialReference> m_poSRS{};
340     bool m_bWriteGDALTags = true;
341     size_t m_nTextLength = 0;
342     mutable std::vector<GUInt64> m_cachedArrayStartIdx{};
343     mutable std::vector<size_t> m_cachedCount{};
344     mutable std::shared_ptr<GDALMDArray> m_poCachedArray{};
345 
346     void ConvertNCToGDAL(GByte*) const;
347     void ConvertGDALToNC(GByte*) const;
348 
349     bool ReadOneElement(const GDALExtendedDataType& src_datatype,
350                         const GDALExtendedDataType& bufferDataType,
351                         const size_t* array_idx,
352                         void* pDstBuffer) const;
353 
354     bool WriteOneElement(const GDALExtendedDataType& dst_datatype,
355                          const GDALExtendedDataType& bufferDataType,
356                          const size_t* array_idx,
357                          const void* pSrcBuffer) const;
358 
359     template< typename BufferType,
360               typename NCGetPutVar1FuncType,
361               typename ReadOrWriteOneElementType >
362     bool IReadWriteGeneric(const size_t* arrayStartIdx,
363                                   const size_t* count,
364                                   const GInt64* arrayStep,
365                                   const GPtrDiff_t* bufferStride,
366                                   const GDALExtendedDataType& bufferDataType,
367                                   BufferType buffer,
368                                   NCGetPutVar1FuncType NCGetPutVar1Func,
369                                   ReadOrWriteOneElementType ReadOrWriteOneElement) const;
370 
371     template< typename BufferType,
372               typename NCGetPutVar1FuncType,
373               typename NCGetPutVaraFuncType,
374               typename NCGetPutVarmFuncType,
375               typename ReadOrWriteOneElementType >
376     bool IReadWrite(const bool bIsRead,
377                     const GUInt64* arrayStartIdx,
378                                     const size_t* count,
379                                     const GInt64* arrayStep,
380                                     const GPtrDiff_t* bufferStride,
381                                     const GDALExtendedDataType& bufferDataType,
382                                     BufferType buffer,
383                                     NCGetPutVar1FuncType NCGetPutVar1Func,
384                                     NCGetPutVaraFuncType NCGetPutVaraFunc,
385                                     NCGetPutVarmFuncType NCGetPutVarmFunc,
386                                     ReadOrWriteOneElementType ReadOrWriteOneElement) const;
387 
388 protected:
389     netCDFVariable(const std::shared_ptr<netCDFSharedResources>& poShared,
390                    int gid, int varid,
391                    const std::vector<std::shared_ptr<GDALDimension>>& dims,
392                    CSLConstList papszOptions);
393 
394     bool IRead(const GUInt64* arrayStartIdx,     // array of size GetDimensionCount()
395                       const size_t* count,                 // array of size GetDimensionCount()
396                       const GInt64* arrayStep,        // step in elements
397                       const GPtrDiff_t* bufferStride, // stride in elements
398                       const GDALExtendedDataType& bufferDataType,
399                       void* pDstBuffer) const override;
400 
401     bool IWrite(const GUInt64* arrayStartIdx,     // array of size GetDimensionCount()
402                       const size_t* count,                 // array of size GetDimensionCount()
403                       const GInt64* arrayStep,        // step in elements
404                       const GPtrDiff_t* bufferStride, // stride in elements
405                       const GDALExtendedDataType& bufferDataType,
406                       const void* pSrcBuffer) override;
407 
408     bool IAdviseRead(const GUInt64* arrayStartIdx,
409                      const size_t* count) const override;
410 
411 public:
Create(const std::shared_ptr<netCDFSharedResources> & poShared,int gid,int varid,const std::vector<std::shared_ptr<GDALDimension>> & dims,CSLConstList papszOptions,bool bCreate)412     static std::shared_ptr<netCDFVariable> Create(
413                    const std::shared_ptr<netCDFSharedResources>& poShared,
414                    int gid, int varid,
415                    const std::vector<std::shared_ptr<GDALDimension>>& dims,
416                    CSLConstList papszOptions,
417                    bool bCreate)
418     {
419         auto var(std::shared_ptr<netCDFVariable>(new netCDFVariable(
420             poShared, gid, varid, dims, papszOptions)));
421         var->SetSelf(var);
422         var->m_bHasWrittenData = !bCreate;
423         return var;
424     }
425 
IsWritable() const426     bool IsWritable() const override { return !m_poShared->IsReadOnly(); }
427 
428     const std::vector<std::shared_ptr<GDALDimension>>& GetDimensions() const override;
429 
430     const GDALExtendedDataType &GetDataType() const override;
431 
432     std::shared_ptr<GDALAttribute> GetAttribute(const std::string& osName) const override;
433 
434     std::vector<std::shared_ptr<GDALAttribute>> GetAttributes(CSLConstList papszOptions) const override;
435 
436     std::shared_ptr<GDALAttribute> CreateAttribute(
437         const std::string& osName,
438         const std::vector<GUInt64>& anDimensions,
439         const GDALExtendedDataType& oDataType,
440         CSLConstList papszOptions) override;
441 
442     const void* GetRawNoDataValue() const override;
443 
444     bool SetRawNoDataValue(const void*) override;
445 
446     std::vector<GUInt64> GetBlockSize() const override;
447 
448     CSLConstList GetStructuralInfo() const override;
449 
GetUnit() const450     const std::string& GetUnit() const override { return m_osUnit; }
451 
452     bool SetUnit(const std::string& osUnit) override;
453 
454     std::shared_ptr<OGRSpatialReference> GetSpatialRef() const override;
455 
456     bool SetSpatialRef(const OGRSpatialReference* poSRS) override;
457 
458     double GetOffset(bool* pbHasOffset, GDALDataType* peStorageType) const override;
459 
460     double GetScale(bool* pbHasScale, GDALDataType* peStorageType) const override;
461 
462     bool SetOffset(double dfOffset, GDALDataType eStorageType) override;
463 
464     bool SetScale(double dfScale, GDALDataType eStorageType) override;
465 
GetGroupId() const466     int GetGroupId() const { return m_gid; }
GetVarId() const467     int GetVarId() const { return m_varid; }
468 
retrieveName(int gid,int varid)469     static std::string retrieveName(int gid, int varid)
470     {
471         CPLMutexHolderD(&hNCMutex);
472         char szName[NC_MAX_NAME + 1] = {};
473         NCDF_ERR(nc_inq_varname(gid, varid, szName));
474         return szName;
475     }
476 };
477 
478 /************************************************************************/
479 /*                       ~netCDFSharedResources()                       */
480 /************************************************************************/
481 
~netCDFSharedResources()482 netCDFSharedResources::~netCDFSharedResources()
483 {
484     CPLMutexHolderD(&hNCMutex);
485 
486     if( m_cdfid > 0 )
487     {
488 #ifdef NCDF_DEBUG
489         CPLDebug("GDAL_netCDF", "calling nc_close( %d)", m_cdfid);
490 #endif
491         int status = nc_close(m_cdfid);
492         NCDF_ERR(status);
493     }
494 
495 #ifdef ENABLE_UFFD
496     if( m_pUffdCtx )
497     {
498         NETCDF_UFFD_UNMAP(m_pUffdCtx);
499     }
500 #endif
501 
502 #ifdef ENABLE_NCDUMP
503     if( m_bFileToDestroyAtClosing )
504         VSIUnlink( m_osFilename );
505 #endif
506 }
507 
508 /************************************************************************/
509 /*                     NCDFGetParentGroupName()                         */
510 /************************************************************************/
511 
NCDFGetParentGroupName(int gid)512 static CPLString NCDFGetParentGroupName(int gid)
513 {
514     int nParentGID = 0;
515     if( nc_inq_grp_parent(gid, &nParentGID) != NC_NOERR )
516         return std::string();
517     return NCDFGetGroupFullName(nParentGID);
518 }
519 
520 /************************************************************************/
521 /*                             netCDFGroup()                            */
522 /************************************************************************/
523 
netCDFGroup(const std::shared_ptr<netCDFSharedResources> & poShared,int gid)524 netCDFGroup::netCDFGroup(const std::shared_ptr<netCDFSharedResources>& poShared, int gid):
525     GDALGroup(NCDFGetParentGroupName(gid), retrieveName(gid)),
526     m_poShared(poShared),
527     m_gid(gid)
528 {
529     if( m_gid == m_poShared->GetCDFId() )
530     {
531         int nFormat = 0;
532         nc_inq_format(m_gid, &nFormat);
533         if( nFormat == NC_FORMAT_CLASSIC )
534         {
535             m_aosStructuralInfo.SetNameValue("NC_FORMAT", "CLASSIC");
536         }
537 #ifdef NC_FORMAT_64BIT_OFFSET
538         else if( nFormat == NC_FORMAT_64BIT_OFFSET )
539         {
540             m_aosStructuralInfo.SetNameValue("NC_FORMAT", "64BIT_OFFSET");
541         }
542 #endif
543 #ifdef NC_FORMAT_CDF5
544         else if( nFormat == NC_FORMAT_CDF5 )
545         {
546             m_aosStructuralInfo.SetNameValue("NC_FORMAT", "CDF5");
547         }
548 #endif
549         else if( nFormat == NC_FORMAT_NETCDF4 )
550         {
551             m_aosStructuralInfo.SetNameValue("NC_FORMAT", "NETCDF4");
552         }
553         else if( nFormat == NC_FORMAT_NETCDF4_CLASSIC )
554         {
555             m_aosStructuralInfo.SetNameValue("NC_FORMAT", "NETCDF4_CLASSIC");
556         }
557     }
558 }
559 
560 /************************************************************************/
561 /*                             CreateGroup()                            */
562 /************************************************************************/
563 
CreateGroup(const std::string & osName,CSLConstList)564 std::shared_ptr<GDALGroup> netCDFGroup::CreateGroup(const std::string& osName,
565                                                  CSLConstList /*papszOptions*/)
566 {
567     if( osName.empty() )
568     {
569         CPLError(CE_Failure, CPLE_NotSupported,
570                  "Empty group name not supported");
571         return nullptr;
572     }
573     CPLMutexHolderD(&hNCMutex);
574     m_poShared->SetDefineMode(true);
575     int nSubGroupId = -1;
576     int ret = nc_def_grp(m_gid, osName.c_str(), &nSubGroupId);
577     NCDF_ERR(ret);
578     if( ret != NC_NOERR )
579         return nullptr;
580     return std::make_shared<netCDFGroup>(m_poShared, nSubGroupId);
581 }
582 
583 /************************************************************************/
584 /*                             CreateDimension()                        */
585 /************************************************************************/
586 
CreateDimension(const std::string & osName,const std::string & osType,const std::string &,GUInt64 nSize,CSLConstList papszOptions)587 std::shared_ptr<GDALDimension> netCDFGroup::CreateDimension(const std::string& osName,
588                                                             const std::string& osType,
589                                                             const std::string&,
590                                                             GUInt64 nSize,
591                                                             CSLConstList papszOptions)
592 {
593     const bool bUnlimited = CPLTestBool(
594         CSLFetchNameValueDef(papszOptions, "UNLIMITED", "FALSE"));
595     if( static_cast<size_t>(nSize) != nSize )
596     {
597         CPLError(CE_Failure, CPLE_AppDefined, "Invalid size");
598         return nullptr;
599     }
600     CPLMutexHolderD(&hNCMutex);
601     m_poShared->SetDefineMode(true);
602     int nDimId = -1;
603     NCDF_ERR(nc_def_dim(m_gid, osName.c_str(), static_cast<size_t>(bUnlimited ? 0 : nSize), &nDimId));
604     if( nDimId < 0 )
605         return nullptr;
606     return std::make_shared<netCDFDimension>(
607         m_poShared, m_gid, nDimId, static_cast<size_t>(nSize), osType);
608 }
609 
610 /************************************************************************/
611 /*                     CreateOrGetComplexDataType()                     */
612 /************************************************************************/
613 
CreateOrGetComplexDataType(int gid,GDALDataType eDT)614 static int CreateOrGetComplexDataType(int gid, GDALDataType eDT)
615 {
616     const char* pszName = "";
617     int nSubTypeId = NC_NAT;
618     switch(eDT)
619     {
620         case GDT_CInt16: pszName = "ComplexInt16"; nSubTypeId = NC_SHORT; break;
621         case GDT_CInt32: pszName = "ComplexInt32"; nSubTypeId = NC_INT; break;
622         case GDT_CFloat32: pszName = "ComplexFloat32"; nSubTypeId = NC_FLOAT; break;
623         case GDT_CFloat64: pszName = "ComplexFloat64"; nSubTypeId = NC_DOUBLE; break;
624         default: CPLAssert(false); break;
625     }
626     int nTypeId = NC_NAT;
627     if( nc_inq_typeid(gid, pszName, &nTypeId) == NC_NOERR )
628     {
629         // We could check that the type definition is really the one we want
630         return nTypeId;
631     }
632     const int nDTSize = GDALGetDataTypeSizeBytes(eDT);
633     NCDF_ERR(nc_def_compound(gid, nDTSize, pszName, &nTypeId));
634     if( nTypeId != NC_NAT )
635     {
636         NCDF_ERR(nc_insert_compound(gid, nTypeId, "real", 0, nSubTypeId));
637         NCDF_ERR(nc_insert_compound(gid, nTypeId, "imag", nDTSize / 2, nSubTypeId));
638     }
639     return nTypeId;
640 }
641 
642 /************************************************************************/
643 /*                    CreateOrGetCompoundDataType()                     */
644 /************************************************************************/
645 
646 static int CreateOrGetType(int gid, const GDALExtendedDataType& oType);
647 
CreateOrGetCompoundDataType(int gid,const GDALExtendedDataType & oType)648 static int CreateOrGetCompoundDataType(int gid, const GDALExtendedDataType& oType)
649 {
650     int nTypeId = NC_NAT;
651     if( nc_inq_typeid(gid, oType.GetName().c_str(), &nTypeId) == NC_NOERR )
652     {
653         // We could check that the type definition is really the one we want
654         return nTypeId;
655     }
656     NCDF_ERR(nc_def_compound(gid, oType.GetSize(), oType.GetName().c_str(), &nTypeId));
657     if( nTypeId != NC_NAT )
658     {
659         for( const auto& comp: oType.GetComponents() )
660         {
661             int nSubTypeId = CreateOrGetType(gid, comp->GetType());
662             if( nSubTypeId == NC_NAT )
663                 return NC_NAT;
664             NCDF_ERR(nc_insert_compound(gid, nTypeId, comp->GetName().c_str(),
665                                         comp->GetOffset(), nSubTypeId));
666         }
667     }
668     return nTypeId;
669 }
670 
671 /************************************************************************/
672 /*                         CreateOrGetType()                            */
673 /************************************************************************/
674 
CreateOrGetType(int gid,const GDALExtendedDataType & oType)675 static int CreateOrGetType(int gid, const GDALExtendedDataType& oType)
676 {
677     int nTypeId = NC_NAT;
678     const auto typeClass = oType.GetClass();
679     if( typeClass == GEDTC_NUMERIC )
680     {
681         switch( oType.GetNumericDataType() )
682         {
683             case GDT_Byte: nTypeId = NC_UBYTE; break;
684             case GDT_UInt16: nTypeId = NC_USHORT; break;
685             case GDT_Int16: nTypeId = NC_SHORT; break;
686             case GDT_UInt32: nTypeId = NC_UINT; break;
687             case GDT_Int32: nTypeId = NC_INT; break;
688             case GDT_Float32: nTypeId = NC_FLOAT; break;
689             case GDT_Float64: nTypeId = NC_DOUBLE; break;
690             case GDT_CInt16:
691             case GDT_CInt32:
692             case GDT_CFloat32:
693             case GDT_CFloat64:
694                 nTypeId = CreateOrGetComplexDataType(gid,
695                                                      oType.GetNumericDataType());
696                 break;
697             default:
698                 break;
699         }
700     }
701     else if( typeClass == GEDTC_STRING )
702     {
703         nTypeId = NC_STRING;
704     }
705     else if( typeClass == GEDTC_COMPOUND )
706     {
707         nTypeId = CreateOrGetCompoundDataType(gid, oType);
708     }
709     return nTypeId;
710 }
711 
712 /************************************************************************/
713 /*                            CreateMDArray()                           */
714 /************************************************************************/
715 
CreateMDArray(const std::string & osName,const std::vector<std::shared_ptr<GDALDimension>> & aoDimensions,const GDALExtendedDataType & oType,CSLConstList papszOptions)716 std::shared_ptr<GDALMDArray> netCDFGroup::CreateMDArray(
717     const std::string& osName,
718     const std::vector<std::shared_ptr<GDALDimension>>& aoDimensions,
719     const GDALExtendedDataType& oType,
720     CSLConstList papszOptions)
721 {
722     if( osName.empty() )
723     {
724         CPLError(CE_Failure, CPLE_NotSupported,
725                  "Empty array name not supported");
726         return nullptr;
727     }
728     CPLMutexHolderD(&hNCMutex);
729     m_poShared->SetDefineMode(true);
730     int nVarId = -1;
731     std::vector<int> anDimIds;
732     std::vector<std::shared_ptr<GDALDimension>> dims;
733     for(const auto& dim: aoDimensions )
734     {
735         int nDimId = -1;
736         auto netCDFDim = std::dynamic_pointer_cast<netCDFDimension>(dim);
737         if( netCDFDim )
738         {
739             nDimId = netCDFDim->GetId();
740         }
741         else
742         {
743             if( nc_inq_dimid(m_gid, dim->GetName().c_str(), &nDimId) == NC_NOERR )
744             {
745                 netCDFDim = std::make_shared<netCDFDimension>(
746                     m_poShared, m_gid, nDimId, 0, dim->GetType());
747                 if( netCDFDim->GetSize() != dim->GetSize() )
748                 {
749                     CPLError(CE_Warning, CPLE_AppDefined,
750                              "Dimension %s already exists, "
751                              "but with a size of " CPL_FRMT_GUIB,
752                              dim->GetName().c_str(),
753                              static_cast<GUIntBig>(netCDFDim->GetSize()));
754                 }
755             }
756             else
757             {
758                 netCDFDim = std::dynamic_pointer_cast<netCDFDimension>(
759                     CreateDimension(dim->GetName(),
760                                     dim->GetType(),
761                                     dim->GetDirection(),
762                                     dim->GetSize(), nullptr));
763                 if( !netCDFDim )
764                     return nullptr;
765                 nDimId = netCDFDim->GetId();
766             }
767         }
768         anDimIds.push_back(nDimId);
769         dims.emplace_back(netCDFDim);
770     }
771     int nTypeId = CreateOrGetType(m_gid, oType);
772     if( nTypeId == NC_NAT )
773     {
774         CPLError(CE_Failure, CPLE_NotSupported,
775                  "Unhandled data type");
776         return nullptr;
777     }
778     const char* pszType = CSLFetchNameValueDef(papszOptions, "NC_TYPE", "");
779     if( (EQUAL(pszType, "") || EQUAL(pszType, "NC_CHAR")) &&
780         dims.size() == 1 &&
781         oType.GetClass() == GEDTC_STRING &&
782         oType.GetMaxStringLength() > 0 )
783     {
784         nTypeId = NC_CHAR;
785         auto dimLength = std::dynamic_pointer_cast<netCDFDimension>(
786                     CreateDimension(aoDimensions[0]->GetName() + "_length",
787                                     std::string(),
788                                     std::string(),
789                                     oType.GetMaxStringLength(), nullptr));
790         if( !dimLength )
791             return nullptr;
792         anDimIds.push_back(dimLength->GetId());
793     }
794     else if( EQUAL(pszType, "NC_BYTE") )
795         nTypeId = NC_BYTE;
796     else if( EQUAL(pszType, "NC_INT64") )
797         nTypeId = NC_INT64;
798     else if( EQUAL(pszType, "NC_UINT64") )
799         nTypeId = NC_UINT64;
800     NCDF_ERR(nc_def_var(m_gid, osName.c_str(), nTypeId,
801                         static_cast<int>(anDimIds.size()),
802                         anDimIds.empty() ? nullptr : anDimIds.data(),
803                         &nVarId));
804     if( nVarId < 0 )
805         return nullptr;
806 
807     const char* pszBlockSize = CSLFetchNameValue(papszOptions, "BLOCKSIZE");
808     if( pszBlockSize &&
809         /* ignore for now BLOCKSIZE for 1-dim string variables created as 2-dim */
810         anDimIds.size() == aoDimensions.size() )
811     {
812         auto aszTokens(CPLStringList(CSLTokenizeString2(pszBlockSize, ",", 0)));
813         if( static_cast<size_t>(aszTokens.size()) != aoDimensions.size() )
814         {
815             CPLError(CE_Failure, CPLE_AppDefined,
816                      "Invalid number of values in BLOCKSIZE");
817             return nullptr;
818         }
819         if( !aoDimensions.empty() )
820         {
821             std::vector<size_t> anChunkSize(aoDimensions.size());
822             for( size_t i = 0; i < anChunkSize.size(); ++i )
823             {
824                 anChunkSize[i] = static_cast<size_t>(CPLAtoGIntBig(aszTokens[i]));
825             }
826             int ret = nc_def_var_chunking(m_gid, nVarId, NC_CHUNKED, &anChunkSize[0]);
827             NCDF_ERR(ret);
828             if( ret != NC_NOERR )
829                 return nullptr;
830         }
831     }
832 
833     const char* pszCompress = CSLFetchNameValue(papszOptions, "COMPRESS");
834     if( pszCompress && EQUAL(pszCompress, "DEFLATE") )
835     {
836         int nZLevel = NCDF_DEFLATE_LEVEL;
837         const char* pszZLevel = CSLFetchNameValue(papszOptions, "ZLEVEL");
838         if( pszZLevel != nullptr )
839         {
840             nZLevel = atoi(pszZLevel);
841             if( !(nZLevel >= 1 && nZLevel <= 9) )
842             {
843                 CPLError(CE_Warning, CPLE_IllegalArg,
844                         "ZLEVEL=%s value not recognised, ignoring.", pszZLevel);
845                 nZLevel = NCDF_DEFLATE_LEVEL;
846             }
847         }
848         int ret = nc_def_var_deflate(m_gid, nVarId,
849                                      TRUE /* shuffle */,
850                                      TRUE /* deflate on */,
851                                      nZLevel);
852         NCDF_ERR(ret);
853         if( ret != NC_NOERR )
854             return nullptr;
855     }
856 
857     const char* pszFilter = CSLFetchNameValue(papszOptions, "FILTER");
858     if( pszFilter )
859     {
860 #ifdef NC_EFILTER
861         const auto aosTokens(CPLStringList(CSLTokenizeString2(pszFilter, ",", 0)));
862         if( !aosTokens.empty() )
863         {
864             const unsigned nFilterId = static_cast<unsigned>(
865                 CPLAtoGIntBig(aosTokens[0]));
866             std::vector<unsigned> anParams;
867             for( int i = 1; i < aosTokens.size(); ++i )
868             {
869                 anParams.push_back(
870                     static_cast<unsigned>(CPLAtoGIntBig(aosTokens[i])));
871             }
872             int ret = nc_def_var_filter(m_gid, nVarId, nFilterId,
873                                         anParams.size(), anParams.data());
874             NCDF_ERR(ret);
875             if( ret != NC_NOERR )
876                 return nullptr;
877         }
878 #else
879         CPLError(CE_Failure, CPLE_NotSupported,
880                  "netCDF 4.6 or later needed for FILTER option");
881         return nullptr;
882 #endif
883     }
884 
885     const bool bChecksum = CPLTestBool(
886         CSLFetchNameValueDef(papszOptions, "CHECKSUM", "FALSE"));
887     if( bChecksum )
888     {
889         int ret = nc_def_var_fletcher32(m_gid, nVarId, TRUE);
890         NCDF_ERR(ret);
891         if( ret != NC_NOERR )
892             return nullptr;
893     }
894 
895     return netCDFVariable::Create(m_poShared, m_gid, nVarId, dims, papszOptions,
896                                   true);
897 }
898 
899 /************************************************************************/
900 /*                          CreateAttribute()                           */
901 /************************************************************************/
902 
CreateAttribute(const std::string & osName,const std::vector<GUInt64> & anDimensions,const GDALExtendedDataType & oDataType,CSLConstList papszOptions)903 std::shared_ptr<GDALAttribute> netCDFGroup::CreateAttribute(
904         const std::string& osName,
905         const std::vector<GUInt64>& anDimensions,
906         const GDALExtendedDataType& oDataType,
907         CSLConstList papszOptions)
908 {
909     return netCDFAttribute::Create(m_poShared, m_gid, NC_GLOBAL,
910                                    osName, anDimensions, oDataType,
911                                    papszOptions);
912 }
913 
914 /************************************************************************/
915 /*                            GetGroupNames()                           */
916 /************************************************************************/
917 
GetGroupNames(CSLConstList papszOptions) const918 std::vector<std::string> netCDFGroup::GetGroupNames(CSLConstList papszOptions) const
919 {
920     CPLMutexHolderD(&hNCMutex);
921     int nSubGroups = 0;
922     NCDF_ERR(nc_inq_grps(m_gid, &nSubGroups, nullptr));
923     if( nSubGroups == 0 )
924     {
925         if( EQUAL(CSLFetchNameValueDef(papszOptions, "GROUP_BY", ""),
926                   "SAME_DIMENSION") )
927         {
928             std::vector<std::string> names;
929             std::set<std::string> oSetDimNames;
930             for( const auto& osArrayName: GetMDArrayNames(nullptr) )
931             {
932                 const auto poArray = OpenMDArray(osArrayName, nullptr);
933                 const auto apoDims = poArray->GetDimensions();
934                 if( apoDims.size() == 1 )
935                 {
936                     const auto osDimName = apoDims[0]->GetName();
937                     if( oSetDimNames.find(osDimName) == oSetDimNames.end() )
938                     {
939                         oSetDimNames.insert(osDimName);
940                         names.emplace_back(osDimName);
941                     }
942                 }
943             }
944             return names;
945         }
946         return {};
947     }
948     std::vector<int> anSubGroupdsIds(nSubGroups);
949     NCDF_ERR(nc_inq_grps(m_gid, nullptr, &anSubGroupdsIds[0]));
950     std::vector<std::string> names;
951     names.reserve(nSubGroups);
952     for( const auto& subgid: anSubGroupdsIds )
953     {
954         char szName[NC_MAX_NAME + 1] = {};
955         NCDF_ERR(nc_inq_grpname(subgid, szName));
956         names.emplace_back(szName);
957     }
958     return names;
959 }
960 
961 /************************************************************************/
962 /*                             OpenGroup()                              */
963 /************************************************************************/
964 
OpenGroup(const std::string & osName,CSLConstList papszOptions) const965 std::shared_ptr<GDALGroup> netCDFGroup::OpenGroup(const std::string& osName,
966                                                   CSLConstList papszOptions) const
967 {
968     CPLMutexHolderD(&hNCMutex);
969     int nSubGroups = 0;
970     // This is weird but nc_inq_grp_ncid() succeeds on a single group file.
971     NCDF_ERR(nc_inq_grps(m_gid, &nSubGroups, nullptr));
972     if( nSubGroups == 0 )
973     {
974         if( EQUAL(CSLFetchNameValueDef(papszOptions, "GROUP_BY", ""),
975             "SAME_DIMENSION") )
976         {
977             const auto oCandidateGroupNames = GetGroupNames(papszOptions);
978             for( const auto& osCandidateGroupName: oCandidateGroupNames )
979             {
980                 if( osCandidateGroupName == osName )
981                 {
982                     auto poThisGroup = std::make_shared<netCDFGroup>(m_poShared, m_gid);
983                     return std::make_shared<netCDFVirtualGroupBySameDimension>(
984                         poThisGroup, osName);
985                 }
986             }
987         }
988         return nullptr;
989     }
990     int nSubGroupId = 0;
991     if( nc_inq_grp_ncid(m_gid, osName.c_str(), &nSubGroupId) != NC_NOERR ||
992         nSubGroupId <= 0 )
993         return nullptr;
994     return std::make_shared<netCDFGroup>(m_poShared, nSubGroupId);
995 }
996 
997 /************************************************************************/
998 /*                         GetMDArrayNames()                            */
999 /************************************************************************/
1000 
GetMDArrayNames(CSLConstList papszOptions) const1001 std::vector<std::string> netCDFGroup::GetMDArrayNames(CSLConstList papszOptions) const
1002 {
1003     CPLMutexHolderD(&hNCMutex);
1004     int nVars = 0;
1005     NCDF_ERR(nc_inq_nvars(m_gid, &nVars));
1006     if( nVars == 0 )
1007         return {};
1008     std::vector<int> anVarIds(nVars);
1009     NCDF_ERR(nc_inq_varids(m_gid, nullptr, &anVarIds[0]));
1010     std::vector<std::string> names;
1011     names.reserve(nVars);
1012     const bool bAll = CPLTestBool(CSLFetchNameValueDef(papszOptions, "SHOW_ALL", "NO"));
1013     const bool bZeroDim = bAll || CPLTestBool(CSLFetchNameValueDef(papszOptions, "SHOW_ZERO_DIM", "NO"));
1014     const bool bCoordinates = bAll || CPLTestBool(CSLFetchNameValueDef(papszOptions, "SHOW_COORDINATES", "YES"));
1015     const bool bBounds = bAll || CPLTestBool(CSLFetchNameValueDef(papszOptions, "SHOW_BOUNDS", "YES"));
1016     const bool bIndexing = bAll || CPLTestBool(CSLFetchNameValueDef(papszOptions, "SHOW_INDEXING", "YES"));
1017     const bool bTime = bAll || CPLTestBool(CSLFetchNameValueDef(papszOptions, "SHOW_TIME", "YES"));
1018     std::set<std::string> ignoreList;
1019     if( !bCoordinates || !bBounds )
1020     {
1021         for( const auto& varid: anVarIds )
1022         {
1023             char **papszTokens = nullptr;
1024             if( !bCoordinates )
1025             {
1026                 char *pszTemp = nullptr;
1027                 if( NCDFGetAttr(m_gid, varid, "coordinates", &pszTemp) == CE_None )
1028                     papszTokens = CSLTokenizeString2(pszTemp, " ", 0);
1029                 CPLFree(pszTemp);
1030             }
1031             if( !bBounds )
1032             {
1033                 char *pszTemp = nullptr;
1034                 if( NCDFGetAttr(m_gid, varid, "bounds", &pszTemp) == CE_None &&
1035                     pszTemp != nullptr && !EQUAL(pszTemp, "") )
1036                     papszTokens = CSLAddString( papszTokens, pszTemp );
1037                 CPLFree(pszTemp);
1038             }
1039             for( char** iter = papszTokens; iter && iter[0]; ++iter )
1040                 ignoreList.insert(*iter);
1041             CSLDestroy(papszTokens);
1042         }
1043     }
1044 
1045     const bool bGroupBySameDimension =
1046         EQUAL(CSLFetchNameValueDef(papszOptions, "GROUP_BY", ""),
1047               "SAME_DIMENSION");
1048     for( const auto& varid: anVarIds )
1049     {
1050         int nVarDims = 0;
1051         NCDF_ERR(nc_inq_varndims(m_gid, varid, &nVarDims));
1052         if( nVarDims == 0 &&! bZeroDim )
1053         {
1054             continue;
1055         }
1056         if( nVarDims == 1 && bGroupBySameDimension )
1057         {
1058             continue;
1059         }
1060 
1061         char szName[NC_MAX_NAME + 1] = {};
1062         NCDF_ERR(nc_inq_varname(m_gid, varid, szName));
1063         if( !bIndexing && nVarDims == 1 )
1064         {
1065             int nDimId = 0;
1066             NCDF_ERR(nc_inq_vardimid(m_gid, varid, &nDimId));
1067             char szDimName[NC_MAX_NAME+1] = {};
1068             NCDF_ERR(nc_inq_dimname(m_gid, nDimId, szDimName));
1069             if( strcmp(szDimName, szName) == 0 )
1070             {
1071                 continue;
1072             }
1073         }
1074 
1075         if( !bTime )
1076         {
1077             char *pszTemp = nullptr;
1078             bool bSkip = false;
1079             if( NCDFGetAttr(m_gid, varid, "standard_name", &pszTemp) == CE_None )
1080             {
1081                 bSkip = pszTemp && strcmp(pszTemp, "time") == 0;
1082             }
1083             CPLFree(pszTemp);
1084             if( bSkip )
1085             {
1086                 continue;
1087             }
1088         }
1089 
1090         if( ignoreList.find(szName) == ignoreList.end() )
1091         {
1092             names.emplace_back(szName);
1093         }
1094     }
1095     return names;
1096 }
1097 
1098 /************************************************************************/
1099 /*                           OpenMDArray()                              */
1100 /************************************************************************/
1101 
OpenMDArray(const std::string & osName,CSLConstList) const1102 std::shared_ptr<GDALMDArray> netCDFGroup::OpenMDArray(
1103                                 const std::string& osName, CSLConstList) const
1104 {
1105     CPLMutexHolderD(&hNCMutex);
1106     int nVarId = 0;
1107     if( nc_inq_varid(m_gid, osName.c_str(), &nVarId) != NC_NOERR )
1108         return nullptr;
1109     return netCDFVariable::Create(m_poShared, m_gid, nVarId,
1110                                   std::vector<std::shared_ptr<GDALDimension>>(),
1111                                   nullptr, false);
1112 }
1113 
1114 /************************************************************************/
1115 /*                         GetDimensions()                              */
1116 /************************************************************************/
1117 
GetDimensions(CSLConstList) const1118 std::vector<std::shared_ptr<GDALDimension>> netCDFGroup::GetDimensions(CSLConstList) const
1119 {
1120     CPLMutexHolderD(&hNCMutex);
1121     int nbDims = 0;
1122     NCDF_ERR(nc_inq_ndims(m_gid, &nbDims));
1123     if( nbDims == 0 )
1124         return {};
1125     std::vector<int> dimids(nbDims);
1126     NCDF_ERR(nc_inq_dimids(m_gid, &nbDims, &dimids[0], FALSE));
1127     std::vector<std::shared_ptr<GDALDimension>> res;
1128     for( int i = 0; i < nbDims; i++ )
1129     {
1130         res.emplace_back(std::make_shared<netCDFDimension>(
1131             m_poShared, m_gid, dimids[i], 0, std::string()));
1132     }
1133     return res;
1134 }
1135 
1136 /************************************************************************/
1137 /*                         GetAttribute()                               */
1138 /************************************************************************/
1139 
GetAttribute(const std::string & osName) const1140 std::shared_ptr<GDALAttribute> netCDFGroup::GetAttribute(const std::string& osName) const
1141 {
1142     CPLMutexHolderD(&hNCMutex);
1143     int nAttId = -1;
1144     if( nc_inq_attid(m_gid, NC_GLOBAL, osName.c_str(), &nAttId) != NC_NOERR )
1145         return nullptr;
1146     return netCDFAttribute::Create(m_poShared, m_gid, NC_GLOBAL, osName);
1147 }
1148 
1149 /************************************************************************/
1150 /*                         GetAttributes()                              */
1151 /************************************************************************/
1152 
GetAttributes(CSLConstList) const1153 std::vector<std::shared_ptr<GDALAttribute>> netCDFGroup::GetAttributes(CSLConstList) const
1154 {
1155     CPLMutexHolderD(&hNCMutex);
1156     std::vector<std::shared_ptr<GDALAttribute>> res;
1157     int nbAttr = 0;
1158     NCDF_ERR(nc_inq_varnatts(m_gid, NC_GLOBAL, &nbAttr));
1159     res.reserve(nbAttr);
1160     for( int i = 0; i < nbAttr; i++ )
1161     {
1162         char szAttrName[NC_MAX_NAME + 1];
1163         szAttrName[0] = 0;
1164         NCDF_ERR(nc_inq_attname(m_gid, NC_GLOBAL, i, szAttrName));
1165         if( !EQUAL(szAttrName, "_NCProperties") )
1166         {
1167             res.emplace_back(netCDFAttribute::Create(
1168                 m_poShared, m_gid, NC_GLOBAL, szAttrName));
1169         }
1170     }
1171     return res;
1172 }
1173 
1174 /************************************************************************/
1175 /*                         GetStructuralInfo()                          */
1176 /************************************************************************/
1177 
GetStructuralInfo() const1178 CSLConstList netCDFGroup::GetStructuralInfo() const
1179 {
1180     return m_aosStructuralInfo.List();
1181 }
1182 
1183 /************************************************************************/
1184 /*                   netCDFVirtualGroupBySameDimension()                */
1185 /************************************************************************/
1186 
netCDFVirtualGroupBySameDimension(const std::shared_ptr<netCDFGroup> & poGroup,const std::string & osDimName)1187 netCDFVirtualGroupBySameDimension::netCDFVirtualGroupBySameDimension(
1188                                     const std::shared_ptr<netCDFGroup>& poGroup,
1189                                     const std::string& osDimName):
1190     GDALGroup(poGroup->GetName(), osDimName),
1191     m_poGroup(poGroup),
1192     m_osDimName(osDimName)
1193 {
1194 }
1195 
1196 /************************************************************************/
1197 /*                         GetMDArrayNames()                            */
1198 /************************************************************************/
1199 
GetMDArrayNames(CSLConstList) const1200 std::vector<std::string> netCDFVirtualGroupBySameDimension::GetMDArrayNames(CSLConstList) const
1201 {
1202     const auto srcNames = m_poGroup->GetMDArrayNames(nullptr);
1203     std::vector<std::string> names;
1204     for( const auto& srcName: srcNames )
1205     {
1206         auto poArray = m_poGroup->OpenMDArray(srcName, nullptr);
1207         if( poArray )
1208         {
1209             const auto apoArrayDims = poArray->GetDimensions();
1210             if( apoArrayDims.size() == 1 &&
1211                 apoArrayDims[0]->GetName() == m_osDimName )
1212             {
1213                 names.emplace_back(srcName);
1214             }
1215         }
1216     }
1217     return names;
1218 }
1219 
1220 /************************************************************************/
1221 /*                           OpenMDArray()                              */
1222 /************************************************************************/
1223 
OpenMDArray(const std::string & osName,CSLConstList papszOptions) const1224 std::shared_ptr<GDALMDArray> netCDFVirtualGroupBySameDimension::OpenMDArray(
1225                                             const std::string& osName,
1226                                             CSLConstList papszOptions) const
1227 {
1228     return m_poGroup->OpenMDArray(osName, papszOptions);
1229 }
1230 
1231 /************************************************************************/
1232 /*                           netCDFDimension()                          */
1233 /************************************************************************/
1234 
netCDFDimension(const std::shared_ptr<netCDFSharedResources> & poShared,int cfid,int dimid,size_t nForcedSize,const std::string & osType)1235 netCDFDimension::netCDFDimension(
1236                 const std::shared_ptr<netCDFSharedResources>& poShared, int cfid,
1237                 int dimid, size_t nForcedSize, const std::string& osType):
1238     GDALDimension(NCDFGetGroupFullName(cfid),
1239                   retrieveName(cfid, dimid),
1240                   osType, // type
1241                   std::string(), // direction
1242                   nForcedSize ? nForcedSize : retrieveSize(cfid, dimid)),
1243     m_poShared(poShared),
1244     m_gid(cfid),
1245     m_dimid(dimid)
1246 {
1247     if( m_osType.empty() && nForcedSize == 0 )
1248     {
1249         auto var = std::dynamic_pointer_cast<netCDFVariable>(
1250             GetIndexingVariable());
1251         if( var )
1252         {
1253             const auto gid = var->GetGroupId();
1254             const auto varid = var->GetVarId();
1255             const auto varname = var->GetName().c_str();
1256             if( NCDFIsVarLongitude(gid, varid, varname) ||
1257                 NCDFIsVarProjectionX(gid, varid, varname) )
1258             {
1259                 m_osType = GDAL_DIM_TYPE_HORIZONTAL_X;
1260                 auto attrPositive = var->GetAttribute(CF_UNITS);
1261                 if( attrPositive )
1262                 {
1263                     const auto val = attrPositive->ReadAsString();
1264                     if( val )
1265                     {
1266                         if( EQUAL(val, CF_DEGREES_EAST) )
1267                         {
1268                             m_osDirection = "EAST";
1269                         }
1270                     }
1271                 }
1272             }
1273             else if( NCDFIsVarLatitude(gid, varid, varname) ||
1274                      NCDFIsVarProjectionY(gid, varid, varname) )
1275             {
1276                 m_osType = GDAL_DIM_TYPE_HORIZONTAL_Y;
1277                 auto attrPositive = var->GetAttribute(CF_UNITS);
1278                 if( attrPositive )
1279                 {
1280                     const auto val = attrPositive->ReadAsString();
1281                     if( val )
1282                     {
1283                         if( EQUAL(val, CF_DEGREES_NORTH) )
1284                         {
1285                             m_osDirection = "NORTH";
1286                         }
1287                     }
1288                 }
1289             }
1290             else if( NCDFIsVarVerticalCoord(gid, varid, varname) )
1291             {
1292                 m_osType = GDAL_DIM_TYPE_VERTICAL;
1293                 auto attrPositive = var->GetAttribute("positive");
1294                 if( attrPositive )
1295                 {
1296                     const auto val = attrPositive->ReadAsString();
1297                     if( val )
1298                     {
1299                         if( EQUAL(val, "up") )
1300                         {
1301                             m_osDirection = "UP";
1302                         }
1303                         else if( EQUAL(val, "down") )
1304                         {
1305                             m_osDirection = "DOWN";
1306                         }
1307                     }
1308                 }
1309             }
1310             else if( NCDFIsVarTimeCoord(gid, varid, varname) )
1311             {
1312                 m_osType = GDAL_DIM_TYPE_TEMPORAL;
1313             }
1314         }
1315     }
1316 }
1317 
1318 /************************************************************************/
1319 /*                         GetIndexingVariable()                        */
1320 /************************************************************************/
1321 
1322 namespace {
1323     struct SetIsInGetIndexingVariable
1324     {
1325         netCDFSharedResources* m_poShared;
1326 
SetIsInGetIndexingVariable__anon1533a8780111::SetIsInGetIndexingVariable1327         explicit SetIsInGetIndexingVariable(netCDFSharedResources* poSharedResources): m_poShared(poSharedResources)
1328         {
1329             m_poShared->SetIsInGetIndexingVariable(true);
1330         }
1331 
~SetIsInGetIndexingVariable__anon1533a8780111::SetIsInGetIndexingVariable1332         ~SetIsInGetIndexingVariable()
1333         {
1334             m_poShared->SetIsInGetIndexingVariable(false);
1335         }
1336     };
1337 }
1338 
GetIndexingVariable() const1339 std::shared_ptr<GDALMDArray> netCDFDimension::GetIndexingVariable() const
1340 {
1341     if( m_poShared->GetIsInIndexingVariable() )
1342         return nullptr;
1343 
1344     SetIsInGetIndexingVariable setterIsInGetIndexingVariable(m_poShared.get());
1345 
1346     CPLMutexHolderD(&hNCMutex);
1347 
1348     // First try to find a variable in this group with the same name as the
1349     // dimension
1350     int nVarId = 0;
1351     if( nc_inq_varid(m_gid, GetName().c_str(), &nVarId) == NC_NOERR )
1352     {
1353         int nDims = 0;
1354         NCDF_ERR(nc_inq_varndims(m_gid, nVarId, &nDims));
1355         int nVarType = NC_NAT;
1356         NCDF_ERR(nc_inq_vartype(m_gid, nVarId, &nVarType));
1357         if( nDims == 1 || (nDims == 2 && nVarType == NC_CHAR) )
1358         {
1359             int anDimIds[2] = {};
1360             NCDF_ERR(nc_inq_vardimid(m_gid, nVarId, anDimIds));
1361             if( anDimIds[0] == m_dimid )
1362             {
1363                 if( nDims == 2 )
1364                 {
1365                     // Check that there is no variable with the same of the
1366                     // second dimension.
1367                     char szExtraDim[NC_MAX_NAME+1] = {};
1368                     NCDF_ERR(nc_inq_dimname(m_gid, anDimIds[1], szExtraDim));
1369                     int nUnused;
1370                     if( nc_inq_varid(m_gid, szExtraDim, &nUnused) == NC_NOERR )
1371                         return nullptr;
1372                 }
1373 
1374                 return netCDFVariable::Create(m_poShared, m_gid, nVarId,
1375                     std::vector<std::shared_ptr<GDALDimension>>(),
1376                     nullptr, false);
1377             }
1378         }
1379     }
1380 
1381     // Otherwise explore the variables in this group to find one that has a
1382     // "coordinates" attribute that references this dimension. If so, let's
1383     // return the variable pointed by the value of "coordinates" as the indexing
1384     // variable. This assumes that there is no other variable that would use
1385     // another variable for the matching dimension of its "coordinates".
1386     netCDFGroup oGroup(m_poShared, m_gid);
1387     const auto arrayNames = oGroup.GetMDArrayNames(nullptr);
1388     std::shared_ptr<GDALMDArray> candidateIndexingVariable;
1389     int nCountCandidateIndexingVariable = 0;
1390     for(const auto& arrayName: arrayNames)
1391     {
1392         const auto poArray = oGroup.OpenMDArray(arrayName, nullptr);
1393         const auto poArrayNC = std::dynamic_pointer_cast<netCDFVariable>(poArray);
1394         if( !poArrayNC )
1395             continue;
1396 
1397         const auto apoArrayDims = poArray->GetDimensions();
1398         if( apoArrayDims.size() == 1 )
1399         {
1400             const auto& poArrayDim = apoArrayDims[0];
1401             const auto poArrayDimNC = std::dynamic_pointer_cast<
1402                                             netCDFDimension>(poArrayDim);
1403             if( poArrayDimNC &&
1404                 poArrayDimNC->m_gid == m_gid &&
1405                 poArrayDimNC->m_dimid == m_dimid )
1406             {
1407                 // If the array doesn't have a coordinates variable, but is a 1D
1408                 // array indexed by our dimension, then use it as the indexing
1409                 // variable, provided it is the only such variable.
1410                 if( nCountCandidateIndexingVariable == 0 )
1411                 {
1412                     candidateIndexingVariable = poArray;
1413                 }
1414                 else
1415                 {
1416                     candidateIndexingVariable.reset();
1417                 }
1418                 nCountCandidateIndexingVariable ++;
1419                 continue;
1420             }
1421         }
1422 
1423         const auto poCoordinates = poArray->GetAttribute("coordinates");
1424         if( !(poCoordinates &&
1425               poCoordinates->GetDataType().GetClass() == GEDTC_STRING) )
1426         {
1427             continue;
1428         }
1429 
1430         // Check that the arrays has as many dimensions as its coordinates attribute
1431         const CPLStringList aosCoordinates(
1432             CSLTokenizeString2(poCoordinates->ReadAsString(), " ", 0));
1433         if( apoArrayDims.size() != static_cast<size_t>(aosCoordinates.size()) )
1434             continue;
1435 
1436         for(size_t i = 0; i < apoArrayDims.size(); ++i)
1437         {
1438             const auto& poArrayDim =  apoArrayDims[i];
1439             const auto poArrayDimNC = std::dynamic_pointer_cast<
1440                                     netCDFDimension>(poArrayDim);
1441 
1442             // Check if the array is indexed by the current dimension
1443             if( !(poArrayDimNC &&
1444                   poArrayDimNC->m_gid == m_gid &&
1445                   poArrayDimNC->m_dimid == m_dimid) )
1446             {
1447                 continue;
1448             }
1449 
1450             // Caution: some datasets have their coordinates variables in the
1451             // same order than dimensions (i.e. from slowest varying to
1452             // fastest varying), while others have the coordinates variables
1453             // in the opposite order.
1454             // Assume same order by default, but if we find the first variable
1455             // to be of longitude/X type, then assume the opposite order.
1456             bool coordinatesInSameOrderThanDimensions = true;
1457             if( aosCoordinates.size() > 1 )
1458             {
1459                 int bFirstGroupId = -1;
1460                 int nFirstVarId = -1;
1461                 if( NCDFResolveVar(poArrayNC->GetGroupId(),
1462                                    aosCoordinates[0],
1463                                    &bFirstGroupId,
1464                                    &nVarId,
1465                                    false) == CE_None &&
1466                    (NCDFIsVarLongitude(bFirstGroupId, nFirstVarId, aosCoordinates[0]) ||
1467                     NCDFIsVarProjectionX(bFirstGroupId, nFirstVarId, aosCoordinates[0])) )
1468                 {
1469                     coordinatesInSameOrderThanDimensions = false;
1470                 }
1471             }
1472 
1473             int nIndexingVarGroupId = -1;
1474             int nIndexingVarId = -1;
1475             const size_t nIdxCoordinate = coordinatesInSameOrderThanDimensions ?
1476                                             i : aosCoordinates.size() - 1 - i;
1477             if( NCDFResolveVar(poArrayNC->GetGroupId(),
1478                                aosCoordinates[nIdxCoordinate],
1479                                &nIndexingVarGroupId,
1480                                &nIndexingVarId,
1481                                false) == CE_None )
1482             {
1483                 return netCDFVariable::Create(m_poShared,
1484                     nIndexingVarGroupId, nIndexingVarId,
1485                     std::vector<std::shared_ptr<GDALDimension>>(),
1486                     nullptr, false);
1487             }
1488         }
1489     }
1490 
1491     return candidateIndexingVariable;
1492 }
1493 
1494 /************************************************************************/
1495 /*                          netCDFVariable()                            */
1496 /************************************************************************/
1497 
netCDFVariable(const std::shared_ptr<netCDFSharedResources> & poShared,int gid,int varid,const std::vector<std::shared_ptr<GDALDimension>> & dims,CSLConstList papszOptions)1498 netCDFVariable::netCDFVariable(const std::shared_ptr<netCDFSharedResources>& poShared,
1499                                int gid, int varid,
1500                                const std::vector<std::shared_ptr<GDALDimension>>& dims,
1501                                CSLConstList papszOptions):
1502     GDALAbstractMDArray(NCDFGetGroupFullName(gid), retrieveName(gid, varid)),
1503     GDALMDArray(NCDFGetGroupFullName(gid), retrieveName(gid, varid)),
1504     m_poShared(poShared),
1505     m_gid(gid),
1506     m_varid(varid),
1507     m_dims(dims)
1508 {
1509     NCDF_ERR(nc_inq_varndims(m_gid, m_varid, &m_nDims));
1510     NCDF_ERR(nc_inq_vartype(m_gid, m_varid, &m_nVarType));
1511     if( m_nDims == 2 && m_nVarType == NC_CHAR )
1512     {
1513         int anDimIds[2] = {};
1514         NCDF_ERR(nc_inq_vardimid(m_gid, m_varid, &anDimIds[0]));
1515 
1516         // Check that there is no variable with the same of the
1517         // second dimension.
1518         char szExtraDim[NC_MAX_NAME+1] = {};
1519         NCDF_ERR(nc_inq_dimname(m_gid, anDimIds[1], szExtraDim));
1520         int nUnused;
1521         if( nc_inq_varid(m_gid, szExtraDim, &nUnused) != NC_NOERR )
1522         {
1523             NCDF_ERR(nc_inq_dimlen(m_gid, anDimIds[1], &m_nTextLength));
1524         }
1525     }
1526 
1527     int nShuffle = 0;
1528     int nDeflate = 0;
1529     int nDeflateLevel = 0;
1530     if( nc_inq_var_deflate(m_gid, m_varid, &nShuffle, &nDeflate, &nDeflateLevel) == NC_NOERR )
1531     {
1532         if( nDeflate )
1533         {
1534             m_aosStructuralInfo.SetNameValue("COMPRESS", "DEFLATE");
1535         }
1536     }
1537     auto unit = netCDFVariable::GetAttribute(CF_UNITS);
1538     if( unit )
1539     {
1540         const char* pszVal = unit->ReadAsString();
1541         if( pszVal )
1542             m_osUnit = pszVal;
1543     }
1544     m_bWriteGDALTags = CPLTestBool(
1545         CSLFetchNameValueDef(papszOptions, "WRITE_GDAL_TAGS", "YES"));
1546 }
1547 
1548 /************************************************************************/
1549 /*                             GetDimensions()                          */
1550 /************************************************************************/
1551 
GetDimensions() const1552 const std::vector<std::shared_ptr<GDALDimension>>& netCDFVariable::GetDimensions() const
1553 {
1554     if( m_nDims == 0 || !m_dims.empty() )
1555         return m_dims;
1556     CPLMutexHolderD(&hNCMutex);
1557     std::vector<int> anDimIds(m_nDims);
1558     NCDF_ERR(nc_inq_vardimid(m_gid, m_varid, &anDimIds[0]));
1559     if( m_nDims == 2 && m_nVarType == NC_CHAR && m_nTextLength > 0 )
1560         anDimIds.resize(1);
1561     m_dims.reserve(m_nDims);
1562     for( const auto& dimid: anDimIds )
1563     {
1564         m_dims.emplace_back(std::make_shared<netCDFDimension>(
1565             m_poShared,
1566             m_poShared->GetBelongingGroupOfDim(m_gid, dimid),
1567             dimid, 0, std::string()));
1568     }
1569     return m_dims;
1570 }
1571 
1572 /************************************************************************/
1573 /*                         GetStructuralInfo()                          */
1574 /************************************************************************/
1575 
GetStructuralInfo() const1576 CSLConstList netCDFVariable::GetStructuralInfo() const
1577 {
1578     return m_aosStructuralInfo.List();
1579 }
1580 
1581 /************************************************************************/
1582 /*                          GetComplexDataType()                        */
1583 /************************************************************************/
1584 
GetComplexDataType(int gid,int nVarType)1585 static GDALDataType GetComplexDataType(int gid, int nVarType)
1586 {
1587     //First enquire and check that the number of fields is 2
1588     size_t nfields = 0;
1589     size_t compoundsize = 0;
1590     char szName[NC_MAX_NAME + 1] = {};
1591     if( nc_inq_compound(gid, nVarType, szName, &compoundsize, &nfields) != NC_NOERR)
1592     {
1593         return GDT_Unknown;
1594     }
1595 
1596     if (nfields != 2 || !STARTS_WITH_CI(szName, "complex"))
1597     {
1598         return GDT_Unknown;
1599     }
1600 
1601     //Now check that that two types are the same in the struct.
1602     nc_type field_type1, field_type2;
1603     int field_dims1, field_dims2;
1604     if ( nc_inq_compound_field(gid, nVarType, 0, nullptr, nullptr, &field_type1, &field_dims1, nullptr) != NC_NOERR)
1605     {
1606         return GDT_Unknown;
1607     }
1608 
1609     if ( nc_inq_compound_field(gid, nVarType, 0, nullptr, nullptr, &field_type2, &field_dims2, nullptr) != NC_NOERR)
1610     {
1611         return GDT_Unknown;
1612     }
1613 
1614     if ((field_type1 != field_type2) || (field_dims1 != field_dims2) || (field_dims1 != 0))
1615     {
1616         return GDT_Unknown;
1617     }
1618 
1619     if (field_type1 == NC_SHORT)
1620     {
1621         return GDT_CInt16;
1622     }
1623     else if (field_type1 == NC_INT)
1624     {
1625         return GDT_CInt32;
1626     }
1627     else if (field_type1 == NC_FLOAT)
1628     {
1629         return GDT_CFloat32;
1630     }
1631     else if (field_type1 == NC_DOUBLE)
1632     {
1633         return GDT_CFloat64;
1634     }
1635 
1636     return GDT_Unknown;
1637 }
1638 
1639 /************************************************************************/
1640 /*                       GetCompoundDataType()                          */
1641 /************************************************************************/
1642 
1643 static bool BuildDataType(int gid, int varid, int nVarType,
1644                           std::unique_ptr<GDALExtendedDataType>& dt,
1645                           bool& bPerfectDataTypeMatch);
1646 
GetCompoundDataType(int gid,int nVarType,std::unique_ptr<GDALExtendedDataType> & dt,bool & bPerfectDataTypeMatch)1647 static bool GetCompoundDataType(int gid, int nVarType,
1648                                 std::unique_ptr<GDALExtendedDataType>& dt,
1649                                 bool& bPerfectDataTypeMatch)
1650 {
1651     size_t nfields = 0;
1652     size_t compoundsize = 0;
1653     char szName[NC_MAX_NAME + 1] = {};
1654     if( nc_inq_compound(gid, nVarType, szName, &compoundsize, &nfields) != NC_NOERR)
1655     {
1656         return false;
1657     }
1658     bPerfectDataTypeMatch = true;
1659     std::vector<std::unique_ptr<GDALEDTComponent>> comps;
1660     for( size_t i = 0; i < nfields; i++ )
1661     {
1662         nc_type field_type = 0;
1663         int field_dims = 0;
1664         size_t field_offset = 0;
1665         char field_name[NC_MAX_NAME + 1] = {};
1666         if ( nc_inq_compound_field(gid, nVarType, static_cast<int>(i),
1667                                    field_name, &field_offset,
1668                                    &field_type, &field_dims, nullptr) != NC_NOERR)
1669         {
1670             return false;
1671         }
1672         if( field_dims != 0 )
1673         {
1674             // We don't support that
1675             return false;
1676         }
1677         std::unique_ptr<GDALExtendedDataType> subDt;
1678         bool bSubPerfectDataTypeMatch = false;
1679         if( !BuildDataType(gid, -1, field_type, subDt, bSubPerfectDataTypeMatch) )
1680         {
1681             return false;
1682         }
1683         if( !bSubPerfectDataTypeMatch )
1684         {
1685             CPLError(CE_Failure, CPLE_NotSupported,
1686                      "Non native GDAL type found in a component of a compound type");
1687             return false;
1688         }
1689         auto comp = std::unique_ptr<GDALEDTComponent>(
1690             new GDALEDTComponent(std::string(field_name), field_offset, *subDt));
1691         comps.emplace_back(std::move(comp));
1692     }
1693     dt.reset(new GDALExtendedDataType(GDALExtendedDataType::Create(
1694         szName, compoundsize, std::move(comps))));
1695 
1696     return dt->GetClass() == GEDTC_COMPOUND;
1697 }
1698 
1699 /************************************************************************/
1700 /*                            BuildDataType()                           */
1701 /************************************************************************/
1702 
BuildDataType(int gid,int varid,int nVarType,std::unique_ptr<GDALExtendedDataType> & dt,bool & bPerfectDataTypeMatch)1703 static bool BuildDataType(int gid, int varid, int nVarType,
1704                           std::unique_ptr<GDALExtendedDataType>& dt,
1705                           bool& bPerfectDataTypeMatch)
1706 {
1707     GDALDataType eDataType = GDT_Unknown;
1708     bPerfectDataTypeMatch = false;
1709     if (NCDFIsUserDefinedType(gid, nVarType))
1710     {
1711         nc_type nBaseType = NC_NAT;
1712         int eClass = 0;
1713         nc_inq_user_type(gid, nVarType, nullptr, nullptr,
1714                          &nBaseType, nullptr, &eClass);
1715         if( eClass == NC_COMPOUND )
1716         {
1717             eDataType = GetComplexDataType(gid, nVarType);
1718             if( eDataType != GDT_Unknown )
1719             {
1720                 bPerfectDataTypeMatch = true;
1721                 dt.reset(new GDALExtendedDataType(GDALExtendedDataType::Create(eDataType)));
1722                 return true;
1723             }
1724             else if( GetCompoundDataType(gid, nVarType, dt, bPerfectDataTypeMatch) )
1725             {
1726                 return true;
1727             }
1728             else
1729             {
1730                 CPLError(CE_Failure, CPLE_NotSupported,
1731                             "Unsupported netCDF compound data type encountered.");
1732                 return false;
1733             }
1734         }
1735         else if( eClass == NC_ENUM )
1736         {
1737             nVarType = nBaseType;
1738         }
1739         else if( eClass == NC_VLEN )
1740         {
1741             CPLError(CE_Failure, CPLE_NotSupported,
1742                      "VLen data type not supported");
1743             return false;
1744         }
1745         else if( eClass == NC_OPAQUE )
1746         {
1747             CPLError(CE_Failure, CPLE_NotSupported,
1748                      "Opaque data type not supported");
1749             return false;
1750         }
1751         else
1752         {
1753             CPLError(CE_Failure, CPLE_NotSupported,
1754                      "Unsupported  netCDF data type encountered.");
1755             return false;
1756         }
1757     }
1758 
1759     if( nVarType == NC_STRING )
1760     {
1761         bPerfectDataTypeMatch = true;
1762         dt.reset(new GDALExtendedDataType(GDALExtendedDataType::CreateString()));
1763         return true;
1764     }
1765     else
1766     {
1767         if( nVarType == NC_BYTE )
1768         {
1769             char *pszTemp = nullptr;
1770             bool bSignedData = true;
1771             if( varid >= 0 && NCDFGetAttr(gid, varid, "_Unsigned", &pszTemp) == CE_None )
1772             {
1773                 if( EQUAL(pszTemp, "true") )
1774                     bSignedData = false;
1775                 else if( EQUAL(pszTemp, "false") )
1776                     bSignedData = true;
1777                 CPLFree(pszTemp);
1778             }
1779             if( !bSignedData )
1780             {
1781                 eDataType = GDT_Byte;
1782                 bPerfectDataTypeMatch = true;
1783             }
1784             else
1785             {
1786                 eDataType = GDT_Int16;
1787             }
1788         }
1789         else if( nVarType == NC_CHAR )
1790         {
1791             // Not sure of this
1792             bPerfectDataTypeMatch = true;
1793             eDataType = GDT_Byte;
1794         }
1795         else if( nVarType == NC_SHORT )
1796         {
1797             bPerfectDataTypeMatch = true;
1798             eDataType = GDT_Int16;
1799         }
1800         else if( nVarType == NC_INT )
1801         {
1802             bPerfectDataTypeMatch = true;
1803             eDataType = GDT_Int32;
1804         }
1805         else if( nVarType == NC_FLOAT )
1806         {
1807             bPerfectDataTypeMatch = true;
1808             eDataType = GDT_Float32;
1809         }
1810         else if( nVarType == NC_DOUBLE )
1811         {
1812             bPerfectDataTypeMatch = true;
1813             eDataType = GDT_Float64;
1814         }
1815         else if( nVarType == NC_UBYTE )
1816         {
1817             bPerfectDataTypeMatch = true;
1818             eDataType = GDT_Byte;
1819         }
1820         else if( nVarType == NC_USHORT )
1821         {
1822             bPerfectDataTypeMatch = true;
1823             eDataType = GDT_UInt16;
1824         }
1825         else if( nVarType == NC_UINT )
1826         {
1827             bPerfectDataTypeMatch = true;
1828             eDataType = GDT_UInt32;
1829         }
1830         else if( nVarType == NC_INT64 )
1831             eDataType = GDT_Float64; // approximation
1832         else if( nVarType == NC_UINT64 )
1833             eDataType = GDT_Float64; // approximation
1834         else
1835         {
1836             CPLError(CE_Failure, CPLE_NotSupported,
1837                         "Unsupported netCDF data type encountered.");
1838             return false;
1839         }
1840     }
1841     dt.reset(new GDALExtendedDataType(GDALExtendedDataType::Create(eDataType)));
1842     return true;
1843 }
1844 
1845 /************************************************************************/
1846 /*                             GetDataType()                            */
1847 /************************************************************************/
1848 
GetDataType() const1849 const GDALExtendedDataType &netCDFVariable::GetDataType() const
1850 {
1851     if( m_dt )
1852         return *m_dt;
1853     CPLMutexHolderD(&hNCMutex);
1854 
1855     if( m_nDims == 2 && m_nVarType == NC_CHAR && m_nTextLength > 0 )
1856     {
1857         m_bPerfectDataTypeMatch = true;
1858         m_dt.reset(new GDALExtendedDataType(
1859             GDALExtendedDataType::CreateString(m_nTextLength)));
1860     }
1861     else
1862     {
1863         m_dt.reset(new GDALExtendedDataType(GDALExtendedDataType::Create(GDT_Unknown)));
1864 
1865         BuildDataType(m_gid, m_varid, m_nVarType, m_dt, m_bPerfectDataTypeMatch);
1866     }
1867     return *m_dt;
1868 }
1869 
1870 /************************************************************************/
1871 /*                              SetUnit()                               */
1872 /************************************************************************/
1873 
SetUnit(const std::string & osUnit)1874 bool netCDFVariable::SetUnit(const std::string& osUnit)
1875 {
1876     if( osUnit.empty() )
1877     {
1878         nc_del_att(m_gid, m_varid, CF_UNITS);
1879         return true;
1880     }
1881     auto poUnits(GetAttribute(CF_UNITS));
1882     if( !poUnits )
1883     {
1884         poUnits = CreateAttribute(CF_UNITS, {},
1885                                   GDALExtendedDataType::CreateString(),
1886                                   nullptr);
1887         if( !poUnits )
1888             return false;
1889     }
1890     return poUnits->Write(osUnit.c_str());
1891 }
1892 
1893 /************************************************************************/
1894 /*                            GetSpatialRef()                           */
1895 /************************************************************************/
1896 
GetSpatialRef() const1897 std::shared_ptr<OGRSpatialReference> netCDFVariable::GetSpatialRef() const
1898 {
1899     if( m_bSRSRead )
1900         return m_poSRS;
1901 
1902     m_bSRSRead = true;
1903     netCDFDataset poDS;
1904     poDS.ReadAttributes(m_gid, m_varid);
1905     int iDimX = 0;
1906     int iDimY = 0;
1907     int iCount = 1;
1908     for( const auto& poDim: GetDimensions() )
1909     {
1910         if( poDim->GetType() == GDAL_DIM_TYPE_HORIZONTAL_X )
1911             iDimX = iCount;
1912         else if( poDim->GetType() == GDAL_DIM_TYPE_HORIZONTAL_Y )
1913             iDimY = iCount;
1914         poDS.papszDimName.AddString(poDim->GetName().c_str());
1915         iCount ++;
1916     }
1917     if( (iDimX == 0 || iDimY == 0) && GetDimensionCount() >= 2 )
1918     {
1919         iDimX = static_cast<int>(GetDimensionCount());
1920         iDimY = iDimX - 1;
1921     }
1922     poDS.SetProjectionFromVar(m_gid, m_varid, true);
1923     auto poSRS = poDS.GetSpatialRef();
1924     if( poSRS )
1925     {
1926         m_poSRS.reset(poSRS->Clone());
1927         if( iDimX > 0 && iDimY > 0 )
1928         {
1929             if( m_poSRS->GetDataAxisToSRSAxisMapping() == std::vector<int>{ 2, 1 } )
1930                 m_poSRS->SetDataAxisToSRSAxisMapping({ iDimY, iDimX });
1931             else
1932                 m_poSRS->SetDataAxisToSRSAxisMapping({ iDimX, iDimY });
1933         }
1934     }
1935 
1936     return m_poSRS;
1937 }
1938 
1939 /************************************************************************/
1940 /*                            SetSpatialRef()                           */
1941 /************************************************************************/
1942 
WriteDimAttr(std::shared_ptr<GDALMDArray> poVar,const char * pszAttrName,const char * pszAttrValue)1943 static void WriteDimAttr(std::shared_ptr<GDALMDArray> poVar,
1944                            const char* pszAttrName,
1945                            const char* pszAttrValue)
1946 {
1947     auto poAttr = poVar->GetAttribute(pszAttrName);
1948     if( poAttr )
1949     {
1950         const char* pszVal = poAttr->ReadAsString();
1951         if( pszVal && !EQUAL(pszVal, pszAttrValue) )
1952         {
1953             CPLError(CE_Warning, CPLE_AppDefined,
1954                         "Variable %s has a %s which is %s and not %s",
1955                         poVar->GetName().c_str(),
1956                         pszAttrName,
1957                         pszVal,
1958                         pszAttrValue);
1959         }
1960     }
1961     else
1962     {
1963         poAttr = poVar->CreateAttribute(
1964             pszAttrName, {},
1965             GDALExtendedDataType::CreateString(),
1966             nullptr);
1967         if( poAttr )
1968             poAttr->Write(pszAttrValue);
1969     }
1970 }
1971 
WriteDimAttrs(std::shared_ptr<GDALDimension> dim,const char * pszStandardName,const char * pszLongName,const char * pszUnits)1972 static void WriteDimAttrs(std::shared_ptr<GDALDimension> dim,
1973                            const char* pszStandardName,
1974                            const char* pszLongName,
1975                            const char* pszUnits)
1976 {
1977     auto poVar = dim->GetIndexingVariable();
1978     if( poVar )
1979     {
1980         WriteDimAttr(poVar, CF_STD_NAME, pszStandardName);
1981         WriteDimAttr(poVar, CF_LNG_NAME, pszLongName);
1982         WriteDimAttr(poVar, CF_UNITS, pszUnits);
1983     }
1984     else
1985     {
1986         CPLError(CE_Warning, CPLE_AppDefined,
1987                  "Dimension %s lacks a indexing variable",
1988                  dim->GetName().c_str());
1989     }
1990 }
1991 
SetSpatialRef(const OGRSpatialReference * poSRS)1992 bool netCDFVariable::SetSpatialRef(const OGRSpatialReference* poSRS)
1993 {
1994     m_bSRSRead = false;
1995     m_poSRS.reset();
1996 
1997     if( poSRS == nullptr )
1998     {
1999         nc_del_att(m_gid, m_varid, CF_GRD_MAPPING);
2000         return true;
2001     }
2002 
2003     char *pszCFProjection = nullptr;
2004     int nSRSVarId = NCDFWriteSRSVariable(m_gid, poSRS, &pszCFProjection,
2005                                          m_bWriteGDALTags);
2006     if( nSRSVarId < 0 || pszCFProjection == nullptr )
2007         return false;
2008 
2009     NCDF_ERR(nc_put_att_text(m_gid, m_varid, CF_GRD_MAPPING,
2010                                 strlen(pszCFProjection),
2011                                 pszCFProjection));
2012     CPLFree(pszCFProjection);
2013 
2014     auto apoDims = GetDimensions();
2015     if( poSRS->IsProjected() )
2016     {
2017         bool bWriteX = false;
2018         bool bWriteY = false;
2019         const char *pszUnits = NCDFGetProjectedCFUnit(poSRS);
2020         for( const auto& poDim: apoDims )
2021         {
2022             const char* pszStandardName = nullptr;
2023             const char* pszLongName = nullptr;
2024             if( poDim->GetType() == GDAL_DIM_TYPE_HORIZONTAL_X ||
2025                 EQUAL(poDim->GetName().c_str(), CF_PROJ_X_VAR_NAME) )
2026             {
2027                 pszStandardName = CF_PROJ_X_COORD;
2028                 pszLongName = CF_PROJ_X_COORD_LONG_NAME;
2029                 bWriteX = true;
2030             }
2031             else if( poDim->GetType() == GDAL_DIM_TYPE_HORIZONTAL_Y ||
2032                 EQUAL(poDim->GetName().c_str(), CF_PROJ_Y_VAR_NAME) )
2033             {
2034                 pszStandardName = CF_PROJ_Y_COORD;
2035                 pszLongName = CF_PROJ_Y_COORD_LONG_NAME;
2036                 bWriteY = true;
2037             }
2038             if( pszStandardName && pszLongName )
2039             {
2040                 WriteDimAttrs(poDim, pszStandardName, pszLongName, pszUnits);
2041             }
2042         }
2043         if( !bWriteX && !bWriteY &&
2044             apoDims.size() >= 2 &&
2045             apoDims[apoDims.size()-2]->GetType().empty() &&
2046             apoDims[apoDims.size()-1]->GetType().empty() &&
2047             apoDims[apoDims.size()-2]->GetIndexingVariable() &&
2048             apoDims[apoDims.size()-1]->GetIndexingVariable() )
2049         {
2050             CPLError(CE_Warning, CPLE_AppDefined,
2051                      "Dimensions of variable %s have no type declared. "
2052                      "Assuming the last one is X, and the preceding one Y",
2053                      GetName().c_str());
2054             WriteDimAttrs(apoDims[apoDims.size()-1],
2055                           CF_PROJ_X_COORD, CF_PROJ_X_COORD_LONG_NAME, pszUnits);
2056             WriteDimAttrs(apoDims[apoDims.size()-2],
2057                           CF_PROJ_Y_COORD, CF_PROJ_Y_COORD_LONG_NAME, pszUnits);
2058         }
2059     }
2060     else if( poSRS->IsGeographic() )
2061     {
2062         bool bWriteX = false;
2063         bool bWriteY = false;
2064         for( const auto& poDim: apoDims )
2065         {
2066             const char* pszStandardName = nullptr;
2067             const char* pszLongName = nullptr;
2068             const char* pszUnits = "";
2069             if( poDim->GetType() == GDAL_DIM_TYPE_HORIZONTAL_X ||
2070                 EQUAL(poDim->GetName().c_str(), CF_LONGITUDE_VAR_NAME) )
2071             {
2072                 pszStandardName = CF_LONGITUDE_STD_NAME;
2073                 pszLongName = CF_LONGITUDE_LNG_NAME;
2074                 pszUnits = CF_DEGREES_EAST;
2075                 bWriteX = true;
2076             }
2077             else if( poDim->GetType() == GDAL_DIM_TYPE_HORIZONTAL_Y ||
2078                 EQUAL(poDim->GetName().c_str(), CF_LATITUDE_VAR_NAME) )
2079             {
2080                 pszStandardName = CF_LATITUDE_STD_NAME;
2081                 pszLongName = CF_LATITUDE_LNG_NAME;
2082                 pszUnits = CF_DEGREES_NORTH;
2083                 bWriteY = true;
2084             }
2085             if( pszStandardName && pszLongName )
2086             {
2087                 WriteDimAttrs(poDim, pszStandardName, pszLongName, pszUnits);
2088             }
2089         }
2090         if( !bWriteX && !bWriteY &&
2091             apoDims.size() >= 2 &&
2092             apoDims[apoDims.size()-2]->GetType().empty() &&
2093             apoDims[apoDims.size()-1]->GetType().empty() &&
2094             apoDims[apoDims.size()-2]->GetIndexingVariable() &&
2095             apoDims[apoDims.size()-1]->GetIndexingVariable() )
2096         {
2097             CPLError(CE_Warning, CPLE_AppDefined,
2098                      "Dimensions of variable %s have no type declared. "
2099                      "Assuming the last one is longitude, "
2100                      "and the preceding one latitude", GetName().c_str());
2101             WriteDimAttrs(apoDims[apoDims.size()-1],
2102                           CF_LONGITUDE_STD_NAME,
2103                           CF_LONGITUDE_LNG_NAME, CF_DEGREES_EAST);
2104             WriteDimAttrs(apoDims[apoDims.size()-2],
2105                           CF_LATITUDE_STD_NAME,
2106                           CF_LATITUDE_LNG_NAME, CF_DEGREES_NORTH);
2107         }
2108     }
2109 
2110     return true;
2111 }
2112 
2113 /************************************************************************/
2114 /*                             GetNCTypeSize()                          */
2115 /************************************************************************/
2116 
GetNCTypeSize(const GDALExtendedDataType & dt,bool bPerfectDataTypeMatch,int nAttType)2117 static size_t GetNCTypeSize(const GDALExtendedDataType& dt,
2118                             bool bPerfectDataTypeMatch,
2119                             int nAttType)
2120 {
2121     auto nElementSize = dt.GetSize();
2122     if( !bPerfectDataTypeMatch )
2123     {
2124         if( nAttType == NC_BYTE )
2125         {
2126             CPLAssert( dt.GetNumericDataType() == GDT_Int16 );
2127             nElementSize = sizeof(signed char);
2128         }
2129         else if( nAttType == NC_INT64 )
2130         {
2131             CPLAssert( dt.GetNumericDataType() == GDT_Float64 );
2132             nElementSize = sizeof(GInt64);
2133         }
2134         else if( nAttType == NC_UINT64 )
2135         {
2136             CPLAssert( dt.GetNumericDataType() == GDT_Float64 );
2137             nElementSize = sizeof(GUInt64);
2138         }
2139         else
2140         {
2141             CPLAssert(false);
2142         }
2143     }
2144     return nElementSize;
2145 }
2146 
2147 /************************************************************************/
2148 /*                   ConvertNCStringsToCPLStrings()                     */
2149 /************************************************************************/
2150 
ConvertNCStringsToCPLStrings(GByte * pBuffer,const GDALExtendedDataType & dt)2151 static void ConvertNCStringsToCPLStrings(GByte* pBuffer,
2152                                          const GDALExtendedDataType& dt)
2153 {
2154     switch( dt.GetClass() )
2155     {
2156         case GEDTC_STRING:
2157         {
2158             char* pszStr;
2159             // cppcheck-suppress pointerSize
2160             memcpy(&pszStr, pBuffer, sizeof(char*));
2161             if( pszStr )
2162             {
2163                 char* pszNewStr = VSIStrdup(pszStr);
2164                 nc_free_string(1, &pszStr);
2165                 // cppcheck-suppress pointerSize
2166                 memcpy(pBuffer, &pszNewStr, sizeof(char*));
2167             }
2168             break;
2169         }
2170 
2171         case GEDTC_NUMERIC:
2172         {
2173             break;
2174         }
2175 
2176         case GEDTC_COMPOUND:
2177         {
2178             const auto& comps = dt.GetComponents();
2179             for( const auto& comp: comps )
2180             {
2181                 ConvertNCStringsToCPLStrings(pBuffer + comp->GetOffset(),
2182                                              comp->GetType());
2183             }
2184             break;
2185         }
2186     }
2187 }
2188 
2189 /************************************************************************/
2190 /*                            FreeNCStrings()                           */
2191 /************************************************************************/
2192 
FreeNCStrings(GByte * pBuffer,const GDALExtendedDataType & dt)2193 static void FreeNCStrings(GByte* pBuffer, const GDALExtendedDataType& dt)
2194 {
2195     switch( dt.GetClass() )
2196     {
2197         case GEDTC_STRING:
2198         {
2199             char* pszStr;
2200             // cppcheck-suppress pointerSize
2201             memcpy(&pszStr, pBuffer, sizeof(char*));
2202             if( pszStr )
2203             {
2204                 nc_free_string(1, &pszStr);
2205             }
2206             break;
2207         }
2208 
2209         case GEDTC_NUMERIC:
2210         {
2211             break;
2212         }
2213 
2214         case GEDTC_COMPOUND:
2215         {
2216             const auto& comps = dt.GetComponents();
2217             for( const auto& comp: comps )
2218             {
2219                 FreeNCStrings(pBuffer + comp->GetOffset(), comp->GetType());
2220             }
2221             break;
2222         }
2223     }
2224 }
2225 
2226 /************************************************************************/
2227 /*                          IReadWriteGeneric()                         */
2228 /************************************************************************/
2229 
2230 namespace {
2231 template<typename T> struct GetGByteType {};
2232 template<> struct GetGByteType<void*>
2233 {
2234     typedef GByte* type;
2235 };
2236 template<> struct GetGByteType<const void*>
2237 {
2238     typedef const GByte* type;
2239 };
2240 }
2241 
2242 template< typename BufferType,
2243           typename NCGetPutVar1FuncType,
2244           typename ReadOrWriteOneElementType >
IReadWriteGeneric(const size_t * arrayStartIdx,const size_t * count,const GInt64 * arrayStep,const GPtrDiff_t * bufferStride,const GDALExtendedDataType & bufferDataType,BufferType buffer,NCGetPutVar1FuncType NCGetPutVar1Func,ReadOrWriteOneElementType ReadOrWriteOneElement) const2245 bool netCDFVariable::IReadWriteGeneric(const size_t* arrayStartIdx,
2246                                   const size_t* count,
2247                                   const GInt64* arrayStep,
2248                                   const GPtrDiff_t* bufferStride,
2249                                   const GDALExtendedDataType& bufferDataType,
2250                                   BufferType buffer,
2251                                   NCGetPutVar1FuncType NCGetPutVar1Func,
2252                                   ReadOrWriteOneElementType ReadOrWriteOneElement) const
2253 {
2254     CPLAssert(m_nDims > 0);
2255     std::vector<size_t> array_idx(m_nDims);
2256     std::vector<size_t> stack_count_iters(m_nDims - 1);
2257     typedef typename GetGByteType<BufferType>::type GBytePtrType;
2258     std::vector<GBytePtrType> stack_ptr(m_nDims);
2259     std::vector<GPtrDiff_t> ptr_inc;
2260     ptr_inc.reserve(m_nDims);
2261     const auto& eArrayEDT = GetDataType();
2262     const bool bSameDT = m_bPerfectDataTypeMatch && eArrayEDT == bufferDataType;
2263     const auto nBufferDTSize = bufferDataType.GetSize();
2264     for( int i = 0; i < m_nDims; i++ )
2265     {
2266         ptr_inc.push_back(bufferStride[i] * nBufferDTSize);
2267     }
2268     const auto nDimsMinus1 = m_nDims - 1;
2269     stack_ptr[0] = static_cast<GBytePtrType>(buffer);
2270 
2271     auto lambdaLastDim = [&](GBytePtrType ptr)
2272     {
2273         array_idx[nDimsMinus1] = arrayStartIdx[nDimsMinus1];
2274         size_t nIters = count[nDimsMinus1];
2275         while(true)
2276         {
2277             if( bSameDT )
2278             {
2279                 int ret = NCGetPutVar1Func(m_gid, m_varid, array_idx.data(), ptr);
2280                 NCDF_ERR(ret);
2281                 if( ret != NC_NOERR )
2282                     return false;
2283             }
2284             else
2285             {
2286                 if( !(this->*ReadOrWriteOneElement)(
2287                                          eArrayEDT, bufferDataType,
2288                                          array_idx.data(), ptr) )
2289                     return false;
2290             }
2291             if( (--nIters) == 0 )
2292                 break;
2293             ptr += ptr_inc[nDimsMinus1];
2294             // CPLUnsanitizedAdd needed as arrayStep[] might be negative, and
2295             // thus automatic conversion from negative to big unsigned might
2296             // occur
2297             array_idx[nDimsMinus1] =
2298                 CPLUnsanitizedAdd<size_t>(
2299                     array_idx[nDimsMinus1],
2300                     static_cast<GPtrDiff_t>(arrayStep[nDimsMinus1]));
2301         }
2302         return true;
2303     };
2304 
2305     if( m_nDims == 1 )
2306     {
2307         return lambdaLastDim(stack_ptr[0]);
2308     }
2309     else if( m_nDims == 2)
2310     {
2311         auto nIters = count[0];
2312         array_idx[0] = arrayStartIdx[0];
2313         while(true)
2314         {
2315             if( !lambdaLastDim(stack_ptr[0]) )
2316                 return false;
2317             if( (--nIters) == 0 )
2318                 break;
2319             stack_ptr[0] += ptr_inc[0];
2320             // CPLUnsanitizedAdd needed as arrayStep[] might be negative, and
2321             // thus automatic conversion from negative to big unsigned might
2322             // occur
2323             array_idx[0] = CPLUnsanitizedAdd<size_t>(
2324                 array_idx[0], static_cast<GPtrDiff_t>(arrayStep[0]));
2325         }
2326     }
2327     else if( m_nDims == 3)
2328     {
2329         stack_count_iters[0] = count[0];
2330         array_idx[0] = arrayStartIdx[0];
2331         while(true)
2332         {
2333             auto nIters = count[1];
2334             array_idx[1] = arrayStartIdx[1];
2335             stack_ptr[1] = stack_ptr[0];
2336             while(true)
2337             {
2338                 if( !lambdaLastDim(stack_ptr[1]) )
2339                     return false;
2340                 if( (--nIters) == 0 )
2341                     break;
2342                 stack_ptr[1] += ptr_inc[1];
2343                 // CPLUnsanitizedAdd needed as arrayStep[] might be negative, and
2344                 // thus automatic conversion from negative to big unsigned might
2345                 // occur
2346                 array_idx[1] = CPLUnsanitizedAdd<size_t>(
2347                     array_idx[1], static_cast<GPtrDiff_t>(arrayStep[1]));
2348             }
2349             if( (--stack_count_iters[0]) == 0 )
2350                 break;
2351             stack_ptr[0] += ptr_inc[0];
2352             // CPLUnsanitizedAdd needed as arrayStep[] might be negative, and
2353             // thus automatic conversion from negative to big unsigned might
2354             // occur
2355             array_idx[0] = CPLUnsanitizedAdd<size_t>(
2356                 array_idx[0], static_cast<GPtrDiff_t>(arrayStep[0]));
2357         }
2358     }
2359     else
2360     {
2361         // Implementation valid for nDims >= 3
2362 
2363         int dimIdx = 0;
2364         // Non-recursive implementation. Hence the gotos
2365         // It might be possible to rewrite this without gotos, but I find they
2366         // make it clearer to understand the recursive nature of the code
2367 lbl_start:
2368         if( dimIdx == nDimsMinus1 - 1 )
2369         {
2370             array_idx[dimIdx] = arrayStartIdx[dimIdx];
2371             auto nIters = count[dimIdx];
2372             while(true)
2373             {
2374                 if( !(lambdaLastDim(stack_ptr[dimIdx])) )
2375                     return false;
2376                 if( (--nIters) == 0 )
2377                     break;
2378                 stack_ptr[dimIdx] += ptr_inc[dimIdx];
2379                 // CPLUnsanitizedAdd needed as arrayStep[] might be negative, and
2380                 // thus automatic conversion from negative to big unsigned might
2381                 // occur
2382                 array_idx[dimIdx] = CPLUnsanitizedAdd<size_t>(
2383                     array_idx[dimIdx], static_cast<GPtrDiff_t>(arrayStep[dimIdx]));
2384             }
2385             // If there was a test if( dimIdx > 0 ), that would be valid for nDims == 2
2386             goto lbl_return_to_caller;
2387         }
2388         else
2389         {
2390             array_idx[dimIdx] = arrayStartIdx[dimIdx];
2391             stack_count_iters[dimIdx] = count[dimIdx];
2392             while(true)
2393             {
2394                 // Simulate a recursive call to the next dimension
2395                 // Implicitly save back count and ptr
2396                 dimIdx ++;
2397                 stack_ptr[dimIdx] = stack_ptr[dimIdx-1];
2398                 goto lbl_start;
2399 lbl_return_to_caller:
2400                 dimIdx --;
2401                 if( (--stack_count_iters[dimIdx]) == 0 )
2402                     break;
2403                 stack_ptr[dimIdx] += ptr_inc[dimIdx];
2404                 // CPLUnsanitizedAdd needed as arrayStep[] might be negative, and
2405                 // thus automatic conversion from negative to big unsigned might
2406                 // occur
2407                 array_idx[dimIdx] = CPLUnsanitizedAdd<size_t>(
2408                     array_idx[dimIdx], static_cast<GPtrDiff_t>(arrayStep[dimIdx]));
2409             }
2410             if( dimIdx > 0 )
2411                 goto lbl_return_to_caller;
2412         }
2413     }
2414 
2415     return true;
2416 }
2417 
2418 /************************************************************************/
2419 /*                          CheckNumericDataType()                      */
2420 /************************************************************************/
2421 
CheckNumericDataType(const GDALExtendedDataType & dt)2422 static bool CheckNumericDataType(const GDALExtendedDataType& dt)
2423 {
2424     const auto klass = dt.GetClass();
2425     if( klass == GEDTC_NUMERIC )
2426         return dt.GetNumericDataType() != GDT_Unknown;
2427     if( klass == GEDTC_STRING )
2428         return false;
2429     CPLAssert( klass == GEDTC_COMPOUND );
2430     const auto& comps = dt.GetComponents();
2431     for( const auto& comp: comps )
2432     {
2433         if( !CheckNumericDataType(comp->GetType()) )
2434             return false;
2435     }
2436     return true;
2437 }
2438 
2439 /************************************************************************/
2440 /*                            IReadWrite()                              */
2441 /************************************************************************/
2442 
2443 template< typename BufferType,
2444           typename NCGetPutVar1FuncType,
2445           typename NCGetPutVaraFuncType,
2446           typename NCGetPutVarmFuncType,
2447           typename ReadOrWriteOneElementType >
IReadWrite(const bool bIsRead,const GUInt64 * arrayStartIdx,const size_t * count,const GInt64 * arrayStep,const GPtrDiff_t * bufferStride,const GDALExtendedDataType & bufferDataType,BufferType buffer,NCGetPutVar1FuncType NCGetPutVar1Func,NCGetPutVaraFuncType NCGetPutVaraFunc,NCGetPutVarmFuncType NCGetPutVarmFunc,ReadOrWriteOneElementType ReadOrWriteOneElement) const2448 bool netCDFVariable::IReadWrite(const bool bIsRead,
2449                                 const GUInt64* arrayStartIdx,
2450                                   const size_t* count,
2451                                   const GInt64* arrayStep,
2452                                   const GPtrDiff_t* bufferStride,
2453                                   const GDALExtendedDataType& bufferDataType,
2454                                   BufferType buffer,
2455                                   NCGetPutVar1FuncType NCGetPutVar1Func,
2456                                   NCGetPutVaraFuncType NCGetPutVaraFunc,
2457                                   NCGetPutVarmFuncType NCGetPutVarmFunc,
2458                                   ReadOrWriteOneElementType ReadOrWriteOneElement) const
2459 {
2460     CPLMutexHolderD(&hNCMutex);
2461     m_poShared->SetDefineMode(false);
2462 
2463     const auto& eDT = GetDataType();
2464     std::vector<size_t> startp;
2465     startp.reserve(m_nDims);
2466     bool bUseSlowPath = !m_bPerfectDataTypeMatch &&
2467         !(bIsRead &&
2468           bufferDataType.GetClass() == GEDTC_NUMERIC &&
2469           eDT.GetClass() == GEDTC_NUMERIC &&
2470           bufferDataType.GetSize() >= eDT.GetSize());
2471     for( int i = 0; i < m_nDims; i++ )
2472     {
2473 #if SIZEOF_VOIDP == 4
2474         if( arrayStartIdx[i] > std::numeric_limits<size_t>::max() )
2475             return false;
2476 #endif
2477         startp.push_back(static_cast<size_t>(arrayStartIdx[i]));
2478 
2479 #if SIZEOF_VOIDP == 4
2480         if( arrayStep[i] < std::numeric_limits<ptrdiff_t>::min() ||
2481             arrayStep[i] > std::numeric_limits<ptrdiff_t>::max() )
2482         {
2483             return false;
2484         }
2485 #endif
2486 
2487         if( count[i] != 1 && arrayStep[i] <= 0 )
2488             bUseSlowPath = true; // netCDF rejects negative or NULL strides
2489 
2490         if( bufferStride[i] < 0 )
2491             bUseSlowPath = true; // and it seems to silently cast to size_t imapp
2492     }
2493 
2494     if( eDT.GetClass() == GEDTC_STRING &&
2495         bufferDataType.GetClass() == GEDTC_STRING &&
2496         m_nVarType == NC_STRING )
2497     {
2498         if( m_nDims == 0 )
2499         {
2500             return (this->*ReadOrWriteOneElement)(eDT, bufferDataType,
2501                                 nullptr, buffer);
2502         }
2503 
2504         return IReadWriteGeneric(
2505                             startp.data(), count, arrayStep,
2506                             bufferStride, bufferDataType, buffer,
2507                             NCGetPutVar1Func, ReadOrWriteOneElement);
2508     }
2509 
2510     if( !CheckNumericDataType(eDT) )
2511         return false;
2512     if( !CheckNumericDataType(bufferDataType) )
2513         return false;
2514 
2515     if( m_nDims == 0 )
2516     {
2517         return (this->*ReadOrWriteOneElement)(eDT, bufferDataType,
2518                               nullptr, buffer);
2519     }
2520 
2521     if( !bUseSlowPath &&
2522         ((GDALDataTypeIsComplex(bufferDataType.GetNumericDataType()) ||
2523          bufferDataType.GetClass() == GEDTC_COMPOUND) &&
2524         bufferDataType == eDT) )
2525     {
2526         // nc_get_varm() not supported for non-atomic types.
2527         ptrdiff_t nExpectedBufferStride = 1;
2528         for( int i = m_nDims; i != 0; )
2529         {
2530             --i;
2531             if( count[i] != 1 &&
2532                 (arrayStep[i] != 1 || bufferStride[i] != nExpectedBufferStride) )
2533             {
2534                 bUseSlowPath = true;
2535                 break;
2536             }
2537             nExpectedBufferStride *= count[i];
2538         }
2539         if( !bUseSlowPath )
2540         {
2541             int ret = NCGetPutVaraFunc(m_gid, m_varid, startp.data(), count, buffer);
2542             NCDF_ERR(ret);
2543             return ret == NC_NOERR;
2544         }
2545     }
2546 
2547     if( bUseSlowPath ||
2548         bufferDataType.GetClass() == GEDTC_COMPOUND ||
2549         eDT.GetClass() == GEDTC_COMPOUND ||
2550         (!bIsRead && bufferDataType.GetNumericDataType() != eDT.GetNumericDataType()) ||
2551         (bIsRead && bufferDataType.GetSize() < eDT.GetSize()) )
2552     {
2553         return IReadWriteGeneric(
2554                             startp.data(), count, arrayStep,
2555                             bufferStride, bufferDataType, buffer,
2556                             NCGetPutVar1Func, ReadOrWriteOneElement);
2557     }
2558 
2559     bUseSlowPath = false;
2560     ptrdiff_t nExpectedBufferStride = 1;
2561     for( int i = m_nDims; i != 0; )
2562     {
2563         --i;
2564         if( count[i] != 1 &&
2565             (arrayStep[i] != 1 || bufferStride[i] != nExpectedBufferStride) )
2566         {
2567             bUseSlowPath = true;
2568             break;
2569         }
2570         nExpectedBufferStride *= count[i];
2571     }
2572     if( !bUseSlowPath )
2573     {
2574         // nc_get_varm() is terribly inefficient, so use nc_get_vara()
2575         // when possible.
2576         int ret = NCGetPutVaraFunc(m_gid, m_varid, startp.data(), count, buffer);
2577         if( ret != NC_NOERR )
2578             return false;
2579         if( bIsRead &&
2580             (!m_bPerfectDataTypeMatch ||
2581              bufferDataType.GetNumericDataType() != eDT.GetNumericDataType()) )
2582         {
2583             // If the buffer data type is "larger" or of the same size as the
2584             // native data type, we can do a in-place conversion
2585             GByte* pabyBuffer = static_cast<GByte*>(const_cast<void*>(buffer));
2586             CPLAssert( bufferDataType.GetSize() >= eDT.GetSize() );
2587             const auto nDTSize = eDT.GetSize();
2588             const auto nBufferDTSize = bufferDataType.GetSize();
2589             if( !m_bPerfectDataTypeMatch && (m_nVarType == NC_CHAR || m_nVarType == NC_BYTE) )
2590             {
2591                 // native NC type translates into GDAL data type of larger size
2592                 for( ptrdiff_t i = nExpectedBufferStride - 1; i >= 0; --i )
2593                 {
2594                     GByte abySrc[2];
2595                     abySrc[0] = *(pabyBuffer + i);
2596                     ConvertNCToGDAL(&abySrc[0]);
2597                     GDALExtendedDataType::CopyValue(
2598                         &abySrc[0], eDT,
2599                         pabyBuffer + i * nBufferDTSize, bufferDataType);
2600                 }
2601             }
2602             else if( !m_bPerfectDataTypeMatch )
2603             {
2604                 // native NC type translates into GDAL data type of same size
2605                 CPLAssert( m_nVarType == NC_INT64 || m_nVarType == NC_UINT64 );
2606                 for( ptrdiff_t i = nExpectedBufferStride - 1; i >= 0; --i )
2607                 {
2608                     ConvertNCToGDAL(pabyBuffer + i * nDTSize);
2609                     GDALExtendedDataType::CopyValue(
2610                         pabyBuffer + i * nDTSize, eDT,
2611                         pabyBuffer + i * nBufferDTSize, bufferDataType);
2612                 }
2613             }
2614             else
2615             {
2616                 for( ptrdiff_t i = nExpectedBufferStride - 1; i >= 0; --i )
2617                 {
2618                     GDALExtendedDataType::CopyValue(
2619                         pabyBuffer + i * nDTSize, eDT,
2620                         pabyBuffer + i * nBufferDTSize, bufferDataType);
2621                 }
2622             }
2623         }
2624         return true;
2625     }
2626     else
2627     {
2628         if( bufferDataType.GetNumericDataType() != eDT.GetNumericDataType() )
2629         {
2630             return IReadWriteGeneric(
2631                                 startp.data(), count, arrayStep,
2632                                 bufferStride, bufferDataType, buffer,
2633                                 NCGetPutVar1Func, ReadOrWriteOneElement);
2634         }
2635         std::vector<ptrdiff_t> stridep;
2636         stridep.reserve(m_nDims);
2637         std::vector<ptrdiff_t> imapp;
2638         imapp.reserve(m_nDims);
2639         for( int i = 0; i < m_nDims; i++ )
2640         {
2641             stridep.push_back(static_cast<ptrdiff_t>(count[i] == 1 ? 1 : arrayStep[i]));
2642             imapp.push_back(static_cast<ptrdiff_t>(bufferStride[i]));
2643         }
2644 
2645         if( !m_poShared->GetImappIsInElements() )
2646         {
2647             const size_t nMul = GetNCTypeSize(eDT,
2648                                             m_bPerfectDataTypeMatch, m_nVarType);
2649             for( int i = 0; i < m_nDims; ++i )
2650             {
2651                 imapp[i] = static_cast<ptrdiff_t>(imapp[i] * nMul);
2652             }
2653         }
2654         int ret = NCGetPutVarmFunc(m_gid, m_varid, startp.data(), count,
2655                             stridep.data(), imapp.data(),
2656                             buffer);
2657         NCDF_ERR(ret);
2658         return ret == NC_NOERR;
2659     }
2660 }
2661 
2662 /************************************************************************/
2663 /*                          ConvertNCToGDAL()                           */
2664 /************************************************************************/
2665 
ConvertNCToGDAL(GByte * buffer) const2666 void netCDFVariable::ConvertNCToGDAL(GByte* buffer) const
2667 {
2668     if( !m_bPerfectDataTypeMatch )
2669     {
2670         if( m_nVarType == NC_CHAR || m_nVarType == NC_BYTE )
2671         {
2672             short s = reinterpret_cast<signed char*>(buffer)[0];
2673             memcpy(buffer, &s, sizeof(s));
2674         }
2675         else if( m_nVarType == NC_INT64 )
2676         {
2677             double v = static_cast<double>(
2678                 reinterpret_cast<GInt64*>(buffer)[0]);
2679             memcpy(buffer, &v, sizeof(v));
2680         }
2681         else if( m_nVarType == NC_UINT64 )
2682         {
2683             double v = static_cast<double>(
2684                 reinterpret_cast<GUInt64*>(buffer)[0]);
2685             memcpy(buffer, &v, sizeof(v));
2686         }
2687     }
2688 }
2689 
2690 /************************************************************************/
2691 /*                           ReadOneElement()                           */
2692 /************************************************************************/
2693 
ReadOneElement(const GDALExtendedDataType & src_datatype,const GDALExtendedDataType & bufferDataType,const size_t * array_idx,void * pDstBuffer) const2694 bool netCDFVariable::ReadOneElement(const GDALExtendedDataType& src_datatype,
2695                                     const GDALExtendedDataType& bufferDataType,
2696                                     const size_t* array_idx,
2697                                     void* pDstBuffer) const
2698 {
2699     if( src_datatype.GetClass() == GEDTC_STRING )
2700     {
2701         char* pszStr = nullptr;
2702         int ret = nc_get_var1_string(m_gid, m_varid, array_idx, &pszStr);
2703         NCDF_ERR(ret);
2704         if( ret != NC_NOERR )
2705             return false;
2706         nc_free_string(1, &pszStr);
2707         GDALExtendedDataType::CopyValue(&pszStr, src_datatype, pDstBuffer, bufferDataType);
2708         return true;
2709     }
2710 
2711     std::vector<GByte> abySrc(std::max(src_datatype.GetSize(),
2712                                 GetNCTypeSize(src_datatype,
2713                                     m_bPerfectDataTypeMatch,
2714                                     m_nVarType)));
2715 
2716     int ret = nc_get_var1(m_gid, m_varid, array_idx, &abySrc[0]);
2717     NCDF_ERR(ret);
2718     if( ret != NC_NOERR )
2719         return false;
2720 
2721     ConvertNCToGDAL(&abySrc[0]);
2722 
2723     GDALExtendedDataType::CopyValue(&abySrc[0], src_datatype, pDstBuffer, bufferDataType);
2724     return true;
2725 }
2726 
2727 /************************************************************************/
2728 /*                                   IRead()                            */
2729 /************************************************************************/
2730 
IRead(const GUInt64 * arrayStartIdx,const size_t * count,const GInt64 * arrayStep,const GPtrDiff_t * bufferStride,const GDALExtendedDataType & bufferDataType,void * pDstBuffer) const2731 bool netCDFVariable::IRead(const GUInt64* arrayStartIdx,
2732                                const size_t* count,
2733                                const GInt64* arrayStep,
2734                                const GPtrDiff_t* bufferStride,
2735                                const GDALExtendedDataType& bufferDataType,
2736                                void* pDstBuffer) const
2737 {
2738     if( m_nDims == 2 && m_nVarType == NC_CHAR && GetDimensions().size() == 1 )
2739     {
2740         CPLMutexHolderD(&hNCMutex);
2741         m_poShared->SetDefineMode(false);
2742 
2743         if( bufferDataType.GetClass() != GEDTC_STRING )
2744             return false;
2745         GByte* pabyDstBuffer = static_cast<GByte*>(pDstBuffer);
2746         size_t array_idx[2] = { static_cast<size_t>(arrayStartIdx[0]), 0 };
2747         size_t array_count[2] = { 1, m_nTextLength };
2748         std::string osTmp(m_nTextLength, 0);
2749         const char* pszTmp = osTmp.c_str();
2750         for(size_t i = 0; i < count[0]; i++ )
2751         {
2752             int ret = nc_get_vara(m_gid, m_varid, array_idx, array_count, &osTmp[0]);
2753             NCDF_ERR(ret);
2754             if( ret != NC_NOERR )
2755                 return false;
2756             // coverity[use_after_free]
2757             GDALExtendedDataType::CopyValue(&pszTmp, GetDataType(), pabyDstBuffer, GetDataType());
2758             array_idx[0] = static_cast<size_t>(array_idx[0] + arrayStep[0]);
2759             pabyDstBuffer += bufferStride[0] * sizeof(char*);
2760         }
2761         return true;
2762     }
2763 
2764     if( m_poCachedArray )
2765     {
2766         const auto nDims = GetDimensionCount();
2767         std::vector<GUInt64> modifiedArrayStartIdx(nDims);
2768         bool canUseCache = true;
2769         for( size_t i = 0; i < nDims; i++ )
2770         {
2771             if( arrayStartIdx[i] >= m_cachedArrayStartIdx[i] &&
2772                 arrayStartIdx[i] + (count[i] - 1) * arrayStep[i] <=
2773                     m_cachedArrayStartIdx[i] + m_cachedCount[i] - 1 )
2774             {
2775                 modifiedArrayStartIdx[i] = arrayStartIdx[i] - m_cachedArrayStartIdx[i];
2776             }
2777             else
2778             {
2779                 canUseCache = false;
2780                 break;
2781             }
2782         }
2783         if( canUseCache )
2784         {
2785             return m_poCachedArray->Read( modifiedArrayStartIdx.data(),
2786                                           count,
2787                                           arrayStep,
2788                                           bufferStride,
2789                                           bufferDataType,
2790                                           pDstBuffer );
2791         }
2792     }
2793 
2794     return IReadWrite
2795                 (true,
2796                  arrayStartIdx, count, arrayStep, bufferStride,
2797                  bufferDataType, pDstBuffer,
2798                  nc_get_var1,
2799                  nc_get_vara,
2800                  nc_get_varm,
2801                  &netCDFVariable::ReadOneElement);
2802 }
2803 
2804 /************************************************************************/
2805 /*                             IAdviseRead()                            */
2806 /************************************************************************/
2807 
IAdviseRead(const GUInt64 * arrayStartIdx,const size_t * count) const2808 bool netCDFVariable::IAdviseRead(const GUInt64* arrayStartIdx,
2809                                  const size_t* count) const
2810 {
2811     const auto nDims = GetDimensionCount();
2812     if( nDims == 0 )
2813         return true;
2814     const auto& eDT = GetDataType();
2815     if( eDT.GetClass() != GEDTC_NUMERIC )
2816         return false;
2817 
2818     auto poMemDriver = static_cast<GDALDriver*>(GDALGetDriverByName("MEM"));
2819     if( poMemDriver == nullptr )
2820         return false;
2821 
2822     m_poCachedArray.reset();
2823 
2824     size_t nElts = 1;
2825     for( size_t i = 0; i < nDims; i++ )
2826         nElts *= count[i];
2827 
2828     void* pData = VSI_MALLOC2_VERBOSE(nElts, eDT.GetSize());
2829     if( pData == nullptr )
2830         return false;
2831 
2832     if( !Read(arrayStartIdx, count, nullptr, nullptr, eDT, pData) )
2833     {
2834         VSIFree(pData);
2835         return false;
2836     }
2837 
2838     auto poDS = poMemDriver->CreateMultiDimensional("", nullptr, nullptr);
2839     auto poGroup = poDS->GetRootGroup();
2840     delete poDS;
2841 
2842     std::vector<std::shared_ptr<GDALDimension>> apoMemDims;
2843     const auto& poDims = GetDimensions();
2844     for( size_t i = 0; i < nDims; i++ )
2845     {
2846         apoMemDims.emplace_back(poGroup->CreateDimension( poDims[i]->GetName(),
2847                                                           std::string(),
2848                                                           std::string(),
2849                                                           count[i],
2850                                                           nullptr ) );
2851     }
2852     m_poCachedArray = poGroup->CreateMDArray(GetName(), apoMemDims, eDT, nullptr);
2853     m_poCachedArray->Write( std::vector<GUInt64>(nDims).data(),
2854                             count,
2855                             nullptr,
2856                             nullptr,
2857                             eDT,
2858                             pData );
2859     m_cachedArrayStartIdx.resize(nDims);
2860     memcpy( &m_cachedArrayStartIdx[0], arrayStartIdx, nDims * sizeof(GUInt64) );
2861     m_cachedCount.resize(nDims);
2862     memcpy( &m_cachedCount[0], count, nDims * sizeof(size_t) );
2863     VSIFree(pData);
2864     return true;
2865 }
2866 
2867 /************************************************************************/
2868 /*                          ConvertGDALToNC()                           */
2869 /************************************************************************/
2870 
ConvertGDALToNC(GByte * buffer) const2871 void netCDFVariable::ConvertGDALToNC(GByte* buffer) const
2872 {
2873     if( !m_bPerfectDataTypeMatch )
2874     {
2875         if( m_nVarType == NC_CHAR || m_nVarType == NC_BYTE )
2876         {
2877             const auto c = static_cast<signed char>(
2878                 reinterpret_cast<short*>(buffer)[0]);
2879             memcpy(buffer, &c, sizeof(c));
2880         }
2881         else if( m_nVarType == NC_INT64 )
2882         {
2883             const auto v = static_cast<GInt64>(
2884                 reinterpret_cast<double*>(buffer)[0]);
2885             memcpy(buffer, &v, sizeof(v));
2886         }
2887         else if( m_nVarType == NC_UINT64 )
2888         {
2889             const auto v = static_cast<GUInt64>(
2890                 reinterpret_cast<double*>(buffer)[0]);
2891             memcpy(buffer, &v, sizeof(v));
2892         }
2893     }
2894 }
2895 
2896 /************************************************************************/
2897 /*                          WriteOneElement()                           */
2898 /************************************************************************/
2899 
WriteOneElement(const GDALExtendedDataType & dst_datatype,const GDALExtendedDataType & bufferDataType,const size_t * array_idx,const void * pSrcBuffer) const2900 bool netCDFVariable::WriteOneElement(const GDALExtendedDataType& dst_datatype,
2901                                     const GDALExtendedDataType& bufferDataType,
2902                                     const size_t* array_idx,
2903                                     const void* pSrcBuffer) const
2904 {
2905     if( dst_datatype.GetClass() == GEDTC_STRING )
2906     {
2907         const char* pszStr = (static_cast<const char* const*>(pSrcBuffer))[0];
2908         int ret = nc_put_var1_string(m_gid, m_varid, array_idx, &pszStr);
2909         NCDF_ERR(ret);
2910         return ret == NC_NOERR;
2911     }
2912 
2913     std::vector<GByte> abyTmp(dst_datatype.GetSize());
2914     GDALExtendedDataType::CopyValue(pSrcBuffer, bufferDataType, &abyTmp[0], dst_datatype);
2915 
2916     ConvertGDALToNC(&abyTmp[0]);
2917 
2918     int ret = nc_put_var1(m_gid, m_varid, array_idx, &abyTmp[0]);
2919     NCDF_ERR(ret);
2920     return ret == NC_NOERR;
2921 }
2922 
2923 /************************************************************************/
2924 /*                                IWrite()                              */
2925 /************************************************************************/
2926 
IWrite(const GUInt64 * arrayStartIdx,const size_t * count,const GInt64 * arrayStep,const GPtrDiff_t * bufferStride,const GDALExtendedDataType & bufferDataType,const void * pSrcBuffer)2927 bool netCDFVariable::IWrite(const GUInt64* arrayStartIdx,
2928                                const size_t* count,
2929                                const GInt64* arrayStep,
2930                                const GPtrDiff_t* bufferStride,
2931                                const GDALExtendedDataType& bufferDataType,
2932                                const void* pSrcBuffer)
2933 {
2934     m_bHasWrittenData = true;
2935 
2936     m_poCachedArray.reset();
2937 
2938     if( m_nDims == 2 && m_nVarType == NC_CHAR && GetDimensions().size() == 1 )
2939     {
2940         CPLMutexHolderD(&hNCMutex);
2941         m_poShared->SetDefineMode(false);
2942 
2943         if( bufferDataType.GetClass() != GEDTC_STRING )
2944             return false;
2945         const char* const* ppszSrcBuffer = static_cast<const char* const*>(pSrcBuffer);
2946         size_t array_idx[2] = { static_cast<size_t>(arrayStartIdx[0]), 0 };
2947         size_t array_count[2] = { 1, m_nTextLength };
2948         std::string osTmp(m_nTextLength, 0);
2949         for(size_t i = 0; i < count[0]; i++ )
2950         {
2951             const char* pszStr = *ppszSrcBuffer;
2952             memset(&osTmp[0], 0, m_nTextLength);
2953             if( pszStr )
2954             {
2955                 size_t nLen = strlen(pszStr);
2956                 memcpy(&osTmp[0], pszStr, std::min(m_nTextLength, nLen));
2957             }
2958             int ret = nc_put_vara(m_gid, m_varid, array_idx, array_count, &osTmp[0]);
2959             NCDF_ERR(ret);
2960             if( ret != NC_NOERR )
2961                 return false;
2962             array_idx[0] = static_cast<size_t>(array_idx[0] + arrayStep[0]);
2963             ppszSrcBuffer += bufferStride[0];
2964         }
2965         return true;
2966     }
2967 
2968     return IReadWrite
2969                 (false,
2970                  arrayStartIdx, count, arrayStep, bufferStride,
2971                  bufferDataType, pSrcBuffer,
2972                  nc_put_var1,
2973                  nc_put_vara,
2974                  nc_put_varm,
2975                  &netCDFVariable::WriteOneElement);
2976 }
2977 
2978 /************************************************************************/
2979 /*                          GetRawNoDataValue()                         */
2980 /************************************************************************/
2981 
GetRawNoDataValue() const2982 const void* netCDFVariable::GetRawNoDataValue() const
2983 {
2984     const auto& dt = GetDataType();
2985     if( m_nVarType == NC_STRING )
2986         return nullptr;
2987 
2988     if( m_bGetRawNoDataValueHasRun )
2989     {
2990         return m_abyNoData.empty() ? nullptr : m_abyNoData.data();
2991     }
2992 
2993     m_bGetRawNoDataValueHasRun = true;
2994     CPLMutexHolderD(&hNCMutex);
2995     std::vector<GByte> abyTmp(std::max(dt.GetSize(),
2996                                 GetNCTypeSize(dt,
2997                                     m_bPerfectDataTypeMatch,
2998                                     m_nVarType)));
2999     int ret = nc_get_att(m_gid, m_varid, _FillValue, &abyTmp[0]);
3000     if( ret != NC_NOERR )
3001     {
3002         m_abyNoData.clear();
3003         return nullptr;
3004     }
3005     ConvertNCToGDAL(&abyTmp[0]);
3006     m_abyNoData.resize(dt.GetSize());
3007     memcpy(&m_abyNoData[0], &abyTmp[0], m_abyNoData.size());
3008     return m_abyNoData.data();
3009 }
3010 
3011 /************************************************************************/
3012 /*                          SetRawNoDataValue()                         */
3013 /************************************************************************/
3014 
SetRawNoDataValue(const void * pNoData)3015 bool netCDFVariable::SetRawNoDataValue(const void* pNoData)
3016 {
3017     GetDataType();
3018     if( m_nVarType == NC_STRING )
3019         return false;
3020 
3021     m_bGetRawNoDataValueHasRun = false;
3022     CPLMutexHolderD(&hNCMutex);
3023     m_poShared->SetDefineMode(true);
3024     int ret;
3025     if( pNoData == nullptr )
3026     {
3027         m_abyNoData.clear();
3028         ret = nc_del_att(m_gid, m_varid, _FillValue);
3029     }
3030     else
3031     {
3032         const auto nSize = GetDataType().GetSize();
3033         m_abyNoData.resize(nSize);
3034         memcpy(&m_abyNoData[0], pNoData, nSize);
3035 
3036         std::vector<GByte> abyTmp(nSize);
3037         memcpy(&abyTmp[0], pNoData, nSize);
3038         ConvertGDALToNC(&abyTmp[0]);
3039 
3040         if( !m_bHasWrittenData )
3041         {
3042             ret = nc_def_var_fill(m_gid, m_varid, NC_FILL, &abyTmp[0]);
3043             NCDF_ERR(ret);
3044         }
3045 
3046         ret = nc_put_att(m_gid, m_varid, _FillValue, m_nVarType, 1, &abyTmp[0]);
3047     }
3048     NCDF_ERR(ret);
3049     if( ret == NC_NOERR )
3050         m_bGetRawNoDataValueHasRun = true;
3051     return ret == NC_NOERR;
3052 }
3053 
3054 /************************************************************************/
3055 /*                               SetScale()                             */
3056 /************************************************************************/
3057 
SetScale(double dfScale,GDALDataType eStorageType)3058 bool netCDFVariable::SetScale(double dfScale, GDALDataType eStorageType)
3059 {
3060     auto poAttr = GetAttribute(CF_SCALE_FACTOR);
3061     if( !poAttr )
3062     {
3063         poAttr = CreateAttribute(CF_SCALE_FACTOR, {},
3064                                  GDALExtendedDataType::Create(
3065                                      eStorageType == GDT_Unknown ? GDT_Float64 : eStorageType),
3066                                  nullptr);
3067     }
3068     if( !poAttr )
3069         return false;
3070     return poAttr->Write(dfScale);
3071 }
3072 
3073 /************************************************************************/
3074 /*                              SetOffset()                             */
3075 /************************************************************************/
3076 
SetOffset(double dfOffset,GDALDataType eStorageType)3077 bool netCDFVariable::SetOffset(double dfOffset, GDALDataType eStorageType)
3078 {
3079     auto poAttr = GetAttribute(CF_ADD_OFFSET);
3080     if( !poAttr )
3081     {
3082         poAttr = CreateAttribute(CF_ADD_OFFSET, {},
3083                                  GDALExtendedDataType::Create(
3084                                      eStorageType == GDT_Unknown ? GDT_Float64 : eStorageType),
3085                                  nullptr);
3086     }
3087     if( !poAttr )
3088         return false;
3089     return poAttr->Write(dfOffset);
3090 }
3091 
3092 /************************************************************************/
3093 /*                               GetScale()                             */
3094 /************************************************************************/
3095 
GetScale(bool * pbHasScale,GDALDataType * peStorageType) const3096 double netCDFVariable::GetScale(bool* pbHasScale, GDALDataType* peStorageType) const
3097 {
3098     auto poAttr = GetAttribute(CF_SCALE_FACTOR);
3099     if( !poAttr || poAttr->GetDataType().GetClass() != GEDTC_NUMERIC )
3100     {
3101         if( pbHasScale )
3102             *pbHasScale = false;
3103         return 1.0;
3104     }
3105     if( pbHasScale )
3106         *pbHasScale = true;
3107     if( peStorageType )
3108         *peStorageType = poAttr->GetDataType().GetNumericDataType();
3109     return poAttr->ReadAsDouble();
3110 }
3111 
3112 /************************************************************************/
3113 /*                               GetOffset()                            */
3114 /************************************************************************/
3115 
GetOffset(bool * pbHasOffset,GDALDataType * peStorageType) const3116 double netCDFVariable::GetOffset(bool* pbHasOffset, GDALDataType* peStorageType) const
3117 {
3118     auto poAttr = GetAttribute(CF_ADD_OFFSET);
3119     if( !poAttr || poAttr->GetDataType().GetClass() != GEDTC_NUMERIC )
3120     {
3121         if( pbHasOffset )
3122             *pbHasOffset = false;
3123         return 0.0;
3124     }
3125     if( pbHasOffset )
3126         *pbHasOffset = true;
3127     if( peStorageType )
3128         *peStorageType = poAttr->GetDataType().GetNumericDataType();
3129     return poAttr->ReadAsDouble();
3130 }
3131 
3132 /************************************************************************/
3133 /*                           GetBlockSize()                             */
3134 /************************************************************************/
3135 
GetBlockSize() const3136 std::vector<GUInt64> netCDFVariable::GetBlockSize() const
3137 {
3138     const auto nDimCount = GetDimensionCount();
3139     std::vector<GUInt64> res(nDimCount);
3140     if( res.empty() )
3141         return res;
3142     int nStorageType = 0;
3143     // We add 1 to the dimension count, for 2D char variables that we
3144     // expose as a 1D variable.
3145     std::vector<size_t> anTemp(1 + nDimCount);
3146     nc_inq_var_chunking(m_gid, m_varid, &nStorageType, &anTemp[0]);
3147     if( nStorageType == NC_CHUNKED )
3148     {
3149         for( size_t i = 0; i < res.size(); ++i )
3150             res[i] = anTemp[i];
3151     }
3152     return res;
3153 }
3154 
3155 /************************************************************************/
3156 /*                         GetAttribute()                               */
3157 /************************************************************************/
3158 
GetAttribute(const std::string & osName) const3159 std::shared_ptr<GDALAttribute> netCDFVariable::GetAttribute(const std::string& osName) const
3160 {
3161     CPLMutexHolderD(&hNCMutex);
3162     int nAttId = -1;
3163     if( nc_inq_attid(m_gid, m_varid, osName.c_str(), &nAttId) != NC_NOERR )
3164         return nullptr;
3165     return netCDFAttribute::Create(m_poShared, m_gid, m_varid, osName);
3166 }
3167 
3168 /************************************************************************/
3169 /*                         GetAttributes()                              */
3170 /************************************************************************/
3171 
GetAttributes(CSLConstList papszOptions) const3172 std::vector<std::shared_ptr<GDALAttribute>> netCDFVariable::GetAttributes(CSLConstList papszOptions) const
3173 {
3174     CPLMutexHolderD(&hNCMutex);
3175     std::vector<std::shared_ptr<GDALAttribute>> res;
3176     int nbAttr = 0;
3177     NCDF_ERR(nc_inq_varnatts(m_gid, m_varid, &nbAttr));
3178     res.reserve(nbAttr);
3179     const bool bShowAll = CPLTestBool(
3180         CSLFetchNameValueDef(papszOptions, "SHOW_ALL", "NO"));
3181     for( int i = 0; i < nbAttr; i++ )
3182     {
3183         char szAttrName[NC_MAX_NAME + 1];
3184         szAttrName[0] = 0;
3185         NCDF_ERR(nc_inq_attname(m_gid, m_varid, i, szAttrName));
3186         if( bShowAll ||
3187             (!EQUAL(szAttrName, _FillValue) &&
3188              !EQUAL(szAttrName, CF_UNITS) &&
3189              !EQUAL(szAttrName, CF_SCALE_FACTOR) &&
3190              !EQUAL(szAttrName, CF_ADD_OFFSET) &&
3191              !EQUAL(szAttrName, CF_GRD_MAPPING) &&
3192              !(EQUAL(szAttrName, "_Unsigned") && m_nVarType == NC_BYTE)) )
3193         {
3194             res.emplace_back(netCDFAttribute::Create(
3195                 m_poShared, m_gid, m_varid, szAttrName));
3196         }
3197     }
3198     return res;
3199 }
3200 
3201 /************************************************************************/
3202 /*                          CreateAttribute()                           */
3203 /************************************************************************/
3204 
CreateAttribute(const std::string & osName,const std::vector<GUInt64> & anDimensions,const GDALExtendedDataType & oDataType,CSLConstList papszOptions)3205 std::shared_ptr<GDALAttribute> netCDFVariable::CreateAttribute(
3206         const std::string& osName,
3207         const std::vector<GUInt64>& anDimensions,
3208         const GDALExtendedDataType& oDataType,
3209         CSLConstList papszOptions)
3210 {
3211     return netCDFAttribute::Create(m_poShared, m_gid, m_varid,
3212                                    osName, anDimensions, oDataType,
3213                                    papszOptions);
3214 }
3215 
3216 /************************************************************************/
3217 /*                    retrieveAttributeParentName()                     */
3218 /************************************************************************/
3219 
retrieveAttributeParentName(int gid,int varid)3220 static CPLString retrieveAttributeParentName(int gid, int varid)
3221 {
3222     auto groupName(NCDFGetGroupFullName(gid));
3223     if( varid == NC_GLOBAL )
3224     {
3225         if( groupName == "/" )
3226             return "/_GLOBAL_";
3227         return groupName + "/_GLOBAL_";
3228     }
3229 
3230     return groupName + "/" +
3231            netCDFVariable::retrieveName(gid, varid);
3232 }
3233 
3234 /************************************************************************/
3235 /*                          netCDFAttribute()                           */
3236 /************************************************************************/
3237 
netCDFAttribute(const std::shared_ptr<netCDFSharedResources> & poShared,int gid,int varid,const std::string & name)3238 netCDFAttribute::netCDFAttribute(
3239                     const std::shared_ptr<netCDFSharedResources>& poShared,
3240                     int gid, int varid, const std::string& name):
3241     GDALAbstractMDArray(retrieveAttributeParentName(gid, varid), name),
3242     GDALAttribute(retrieveAttributeParentName(gid, varid), name),
3243     m_poShared(poShared),
3244     m_gid(gid),
3245     m_varid(varid)
3246 {
3247     CPLMutexHolderD(&hNCMutex);
3248     size_t nLen = 0;
3249     NCDF_ERR(nc_inq_atttype(m_gid, m_varid, GetName().c_str(), &m_nAttType));
3250     NCDF_ERR(nc_inq_attlen(m_gid, m_varid, GetName().c_str(), &nLen));
3251     if( m_nAttType == NC_CHAR )
3252     {
3253         m_nTextLength = nLen;
3254     }
3255     else if( nLen > 1 )
3256     {
3257         m_dims.emplace_back(std::make_shared<GDALDimension>(
3258             std::string(), "length", std::string(), std::string(), nLen));
3259     }
3260 }
3261 
3262 /************************************************************************/
3263 /*                          netCDFAttribute()                           */
3264 /************************************************************************/
3265 
netCDFAttribute(const std::shared_ptr<netCDFSharedResources> & poShared,int gid,int varid,const std::string & osName,const std::vector<GUInt64> & anDimensions,const GDALExtendedDataType & oDataType,CSLConstList papszOptions)3266 netCDFAttribute::netCDFAttribute(const std::shared_ptr<netCDFSharedResources>& poShared,
3267                    int gid, int varid, const std::string& osName,
3268                    const std::vector<GUInt64>& anDimensions,
3269                    const GDALExtendedDataType& oDataType,
3270                    CSLConstList papszOptions):
3271     GDALAbstractMDArray(retrieveAttributeParentName(gid, varid), osName),
3272     GDALAttribute(retrieveAttributeParentName(gid, varid), osName),
3273     m_poShared(poShared),
3274     m_gid(gid),
3275     m_varid(varid)
3276 {
3277     CPLMutexHolderD(&hNCMutex);
3278     m_bPerfectDataTypeMatch = true;
3279     m_nAttType = CreateOrGetType(gid, oDataType);
3280     m_dt.reset(new GDALExtendedDataType(oDataType));
3281     if( !anDimensions.empty() )
3282     {
3283         m_dims.emplace_back(std::make_shared<GDALDimension>(
3284             std::string(), "length", std::string(), std::string(), anDimensions[0]));
3285     }
3286 
3287     const char* pszType = CSLFetchNameValueDef(papszOptions, "NC_TYPE", "");
3288     if( oDataType.GetClass() == GEDTC_STRING && anDimensions.empty() &&
3289         (EQUAL(pszType, "") || EQUAL(pszType, "NC_CHAR")) )
3290     {
3291         m_nAttType = NC_CHAR;
3292     }
3293     else if( oDataType.GetNumericDataType() == GDT_Int16 &&
3294              EQUAL(CSLFetchNameValueDef(papszOptions, "NC_TYPE", ""), "NC_BYTE") )
3295     {
3296         m_bPerfectDataTypeMatch = false;
3297         m_nAttType = NC_BYTE;
3298     }
3299     else if( oDataType.GetNumericDataType() == GDT_Float64 )
3300     {
3301         if( EQUAL(pszType, "NC_INT64") )
3302         {
3303             m_bPerfectDataTypeMatch = false;
3304             m_nAttType = NC_INT64;
3305         }
3306         else if( EQUAL(pszType, "NC_UINT64") )
3307         {
3308             m_bPerfectDataTypeMatch = false;
3309             m_nAttType = NC_UINT64;
3310         }
3311     }
3312 }
3313 
3314 /************************************************************************/
3315 /*                              Create()                                */
3316 /************************************************************************/
3317 
Create(const std::shared_ptr<netCDFSharedResources> & poShared,int gid,int varid,const std::string & name)3318 std::shared_ptr<netCDFAttribute> netCDFAttribute::Create(
3319                     const std::shared_ptr<netCDFSharedResources>& poShared,
3320                     int gid, int varid, const std::string& name)
3321 {
3322     auto attr(std::shared_ptr<netCDFAttribute>(
3323         new netCDFAttribute(poShared, gid, varid, name)));
3324     attr->SetSelf(attr);
3325     return attr;
3326 }
3327 
Create(const std::shared_ptr<netCDFSharedResources> & poShared,int gid,int varid,const std::string & osName,const std::vector<GUInt64> & anDimensions,const GDALExtendedDataType & oDataType,CSLConstList papszOptions)3328 std::shared_ptr<netCDFAttribute> netCDFAttribute::Create(
3329                    const std::shared_ptr<netCDFSharedResources>& poShared,
3330                    int gid, int varid, const std::string& osName,
3331                    const std::vector<GUInt64>& anDimensions,
3332                    const GDALExtendedDataType& oDataType,
3333                    CSLConstList papszOptions)
3334 {
3335     if( poShared->IsReadOnly() )
3336     {
3337         CPLError(CE_Failure, CPLE_AppDefined,
3338                  "CreateAttribute() not supported on read-only file");
3339         return nullptr;
3340     }
3341     if( anDimensions.size() > 1 )
3342     {
3343         CPLError(CE_Failure, CPLE_NotSupported,
3344                  "Only 0 or 1-dimensional attribute are supported");
3345         return nullptr;
3346     }
3347     auto attr(std::shared_ptr<netCDFAttribute>(
3348         new netCDFAttribute(poShared, gid, varid, osName,
3349                             anDimensions, oDataType,
3350                             papszOptions)));
3351     if( attr->m_nAttType == NC_NAT )
3352         return nullptr;
3353     attr->SetSelf(attr);
3354     return attr;
3355 }
3356 
3357 /************************************************************************/
3358 /*                             GetDataType()                            */
3359 /************************************************************************/
3360 
GetDataType() const3361 const GDALExtendedDataType &netCDFAttribute::GetDataType() const
3362 {
3363     if( m_dt )
3364         return *m_dt;
3365     CPLMutexHolderD(&hNCMutex);
3366 
3367     if( m_nAttType == NC_CHAR )
3368     {
3369         m_dt.reset(new GDALExtendedDataType(GDALExtendedDataType::CreateString()));
3370     }
3371     else
3372     {
3373         m_dt.reset(new GDALExtendedDataType(GDALExtendedDataType::Create(GDT_Unknown)));
3374         BuildDataType(m_gid, m_varid, m_nAttType, m_dt, m_bPerfectDataTypeMatch);
3375     }
3376 
3377     return *m_dt;
3378 }
3379 
3380 /************************************************************************/
3381 /*                                   IRead()                            */
3382 /************************************************************************/
3383 
IRead(const GUInt64 * arrayStartIdx,const size_t * count,const GInt64 * arrayStep,const GPtrDiff_t * bufferStride,const GDALExtendedDataType & bufferDataType,void * pDstBuffer) const3384 bool netCDFAttribute::IRead(const GUInt64* arrayStartIdx,
3385                                const size_t* count,
3386                                const GInt64* arrayStep,
3387                                const GPtrDiff_t* bufferStride,
3388                                const GDALExtendedDataType& bufferDataType,
3389                                void* pDstBuffer) const
3390 {
3391     CPLMutexHolderD(&hNCMutex);
3392 
3393     if( m_nAttType == NC_STRING )
3394     {
3395         CPLAssert( GetDataType().GetClass() == GEDTC_STRING );
3396         std::vector<char*> apszStrings(static_cast<size_t>(GetTotalElementsCount()));
3397         int ret = nc_get_att_string(m_gid, m_varid, GetName().c_str(), &apszStrings[0]);
3398         NCDF_ERR(ret);
3399         if( ret != NC_NOERR )
3400             return false;
3401         if( m_dims.empty() )
3402         {
3403             const char* pszStr = apszStrings[0];
3404             GDALExtendedDataType::CopyValue(&pszStr, GetDataType(), pDstBuffer, bufferDataType);
3405         }
3406         else
3407         {
3408             GByte* pabyDstBuffer = static_cast<GByte*>(pDstBuffer);
3409             for(size_t i = 0; i < count[0]; i++ )
3410             {
3411                 auto srcIdx = static_cast<size_t>(arrayStartIdx[0] + arrayStep[0] * i);
3412                 const char* pszStr = apszStrings[srcIdx];
3413                 GDALExtendedDataType::CopyValue(&pszStr, GetDataType(),
3414                             pabyDstBuffer, bufferDataType);
3415                 pabyDstBuffer += sizeof(char*) * bufferStride[0];
3416             }
3417         }
3418         nc_free_string(apszStrings.size(), &apszStrings[0]);
3419         return true;
3420     }
3421 
3422     if( m_nAttType == NC_CHAR )
3423     {
3424         CPLAssert( GetDataType().GetClass() == GEDTC_STRING );
3425         CPLAssert( m_dims.empty() );
3426         if( bufferDataType != GetDataType() )
3427         {
3428             std::string osStr;
3429             osStr.resize(m_nTextLength);
3430             int ret = nc_get_att_text(m_gid, m_varid, GetName().c_str(), &osStr[0]);
3431             NCDF_ERR(ret);
3432             if( ret != NC_NOERR )
3433                 return false;
3434             const char* pszStr = osStr.c_str();
3435             GDALExtendedDataType::CopyValue(&pszStr, GetDataType(), pDstBuffer, bufferDataType);
3436         }
3437         else
3438         {
3439             char* pszStr = static_cast<char*>(CPLCalloc(1, m_nTextLength + 1));
3440             int ret = nc_get_att_text(m_gid, m_varid, GetName().c_str(), pszStr);
3441             NCDF_ERR(ret);
3442             if( ret != NC_NOERR )
3443             {
3444                 CPLFree(pszStr);
3445                 return false;
3446             }
3447             *static_cast<char**>(pDstBuffer) = pszStr;
3448         }
3449         return true;
3450     }
3451 
3452     const auto dt(GetDataType());
3453     if( dt.GetClass() == GEDTC_NUMERIC &&
3454         dt.GetNumericDataType() == GDT_Unknown )
3455     {
3456         return false;
3457     }
3458 
3459     CPLAssert( dt.GetClass() != GEDTC_STRING );
3460     const bool bFastPath =
3461         ((m_dims.size() == 1 &&
3462           arrayStartIdx[0] == 0 &&
3463           count[0] == m_dims[0]->GetSize() &&
3464           arrayStep[0] == 1 &&
3465           bufferStride[0] == 1 ) ||
3466          m_dims.empty() ) &&
3467         m_bPerfectDataTypeMatch &&
3468         bufferDataType == dt &&
3469         dt.GetSize() > 0;
3470     if( bFastPath )
3471     {
3472         int ret = nc_get_att(m_gid, m_varid, GetName().c_str(), pDstBuffer);
3473         NCDF_ERR(ret);
3474         if( ret == NC_NOERR )
3475         {
3476             ConvertNCStringsToCPLStrings(static_cast<GByte*>(pDstBuffer), dt);
3477         }
3478         return ret == NC_NOERR;
3479     }
3480 
3481     const auto nElementSize = GetNCTypeSize(dt,
3482                                             m_bPerfectDataTypeMatch,
3483                                             m_nAttType);
3484     if( nElementSize == 0 )
3485         return false;
3486     const auto nOutputDTSize = bufferDataType.GetSize();
3487     std::vector<GByte> abyBuffer(
3488         static_cast<size_t>(GetTotalElementsCount()) * nElementSize);
3489     int ret = nc_get_att(m_gid, m_varid, GetName().c_str(), &abyBuffer[0]);
3490     NCDF_ERR(ret);
3491     if( ret != NC_NOERR )
3492         return false;
3493 
3494     GByte* pabySrcBuffer = m_dims.empty() ? abyBuffer.data() :
3495         abyBuffer.data() +
3496             static_cast<size_t>(arrayStartIdx[0]) * nElementSize;
3497     GByte* pabyDstBuffer = static_cast<GByte*>(pDstBuffer);
3498     for(size_t i = 0; i < (m_dims.empty() ? 1 : count[0]); i++ )
3499     {
3500         GByte abyTmpBuffer[sizeof(double)];
3501         const GByte* pabySrcElement = pabySrcBuffer;
3502         if( !m_bPerfectDataTypeMatch )
3503         {
3504             if( m_nAttType == NC_BYTE )
3505             {
3506                 short s = reinterpret_cast<const signed char*>(pabySrcBuffer)[0];
3507                 memcpy(abyTmpBuffer, &s, sizeof(s));
3508                 pabySrcElement = abyTmpBuffer;
3509             }
3510             else if( m_nAttType == NC_INT64 )
3511             {
3512                 double v = static_cast<double>(
3513                     reinterpret_cast<const GInt64*>(pabySrcBuffer)[0]);
3514                 memcpy(abyTmpBuffer, &v, sizeof(v));
3515                 pabySrcElement = abyTmpBuffer;
3516             }
3517             else if( m_nAttType == NC_UINT64 )
3518             {
3519                 double v = static_cast<double>(
3520                     reinterpret_cast<const GUInt64*>(pabySrcBuffer)[0]);
3521                 memcpy(abyTmpBuffer, &v, sizeof(v));
3522                 pabySrcElement = abyTmpBuffer;
3523             }
3524             else
3525             {
3526                 CPLAssert(false);
3527             }
3528         }
3529         GDALExtendedDataType::CopyValue(pabySrcElement, dt,
3530                     pabyDstBuffer, bufferDataType);
3531         FreeNCStrings(pabySrcBuffer, dt);
3532         if( !m_dims.empty() )
3533         {
3534             pabySrcBuffer += static_cast<std::ptrdiff_t>(arrayStep[0] * nElementSize);
3535             pabyDstBuffer += nOutputDTSize * bufferStride[0];
3536         }
3537     }
3538 
3539     return true;
3540 }
3541 
3542 /************************************************************************/
3543 /*                                IWrite()                              */
3544 /************************************************************************/
3545 
IWrite(const GUInt64 * arrayStartIdx,const size_t * count,const GInt64 * arrayStep,const GPtrDiff_t * bufferStride,const GDALExtendedDataType & bufferDataType,const void * pSrcBuffer)3546 bool netCDFAttribute::IWrite(const GUInt64* arrayStartIdx,
3547                                const size_t* count,
3548                                const GInt64* arrayStep,
3549                                const GPtrDiff_t* bufferStride,
3550                                const GDALExtendedDataType& bufferDataType,
3551                                const void* pSrcBuffer)
3552 {
3553     CPLMutexHolderD(&hNCMutex);
3554 
3555     if( m_dims.size() == 1 &&
3556         (arrayStartIdx[0] != 0 ||
3557          count[0] != m_dims[0]->GetSize() || arrayStep[0] != 1) )
3558     {
3559         CPLError(CE_Failure, CPLE_NotSupported,
3560                     "Only contiguous writing of attribute values supported");
3561         return false;
3562     }
3563 
3564     m_poShared->SetDefineMode(true);
3565 
3566     if( m_nAttType == NC_STRING )
3567     {
3568         CPLAssert( GetDataType().GetClass() == GEDTC_STRING );
3569         if( m_dims.empty() )
3570         {
3571             char* pszStr = nullptr;
3572             const char* pszStrConst;
3573             if( bufferDataType != GetDataType() )
3574             {
3575                 GDALExtendedDataType::CopyValue(pSrcBuffer, bufferDataType, &pszStr, GetDataType());
3576                 pszStrConst = pszStr;
3577             }
3578             else
3579             {
3580                 memcpy(&pszStrConst, pSrcBuffer, sizeof(const char*));
3581             }
3582             int ret = nc_put_att_string(m_gid, m_varid, GetName().c_str(),
3583                                         1, &pszStrConst);
3584             CPLFree(pszStr);
3585             NCDF_ERR(ret);
3586             if( ret != NC_NOERR )
3587                 return false;
3588             return true;
3589         }
3590 
3591         int ret;
3592         if( bufferDataType != GetDataType() )
3593         {
3594             std::vector<char*> apszStrings(count[0]);
3595             const char** ppszStr;
3596             memcpy(&ppszStr, &pSrcBuffer, sizeof(const char**));
3597             for( size_t i = 0; i < count[0]; i++ )
3598             {
3599                 GDALExtendedDataType::CopyValue(&ppszStr[i], bufferDataType, &apszStrings[i], GetDataType());
3600             }
3601             ret = nc_put_att_string(m_gid, m_varid, GetName().c_str(),
3602                                         count[0],
3603                                         const_cast<const char**>(&apszStrings[0]));
3604             for( size_t i = 0; i < count[0]; i++ )
3605             {
3606                 CPLFree(apszStrings[i]);
3607             }
3608         }
3609         else
3610         {
3611             const char** ppszStr;
3612             memcpy(&ppszStr, &pSrcBuffer, sizeof(const char**));
3613             ret = nc_put_att_string(m_gid, m_varid, GetName().c_str(),
3614                                     count[0],
3615                                     ppszStr);
3616         }
3617         NCDF_ERR(ret);
3618         if( ret != NC_NOERR )
3619             return false;
3620         return true;
3621     }
3622 
3623     if( m_nAttType == NC_CHAR )
3624     {
3625         CPLAssert( GetDataType().GetClass() == GEDTC_STRING );
3626         CPLAssert( m_dims.empty() );
3627         char* pszStr = nullptr;
3628         const char* pszStrConst;
3629         if( bufferDataType != GetDataType() )
3630         {
3631             GDALExtendedDataType::CopyValue(pSrcBuffer, bufferDataType, &pszStr, GetDataType());
3632             pszStrConst = pszStr;
3633         }
3634         else
3635         {
3636             memcpy(&pszStrConst, pSrcBuffer, sizeof(const char*));
3637         }
3638         m_nTextLength = pszStrConst ? strlen(pszStrConst) : 0;
3639         int ret = nc_put_att_text(m_gid, m_varid, GetName().c_str(),
3640                                   m_nTextLength, pszStrConst);
3641         CPLFree(pszStr);
3642         NCDF_ERR(ret);
3643         if( ret != NC_NOERR )
3644             return false;
3645         return true;
3646     }
3647 
3648     const auto dt(GetDataType());
3649     if( dt.GetClass() == GEDTC_NUMERIC &&
3650         dt.GetNumericDataType() == GDT_Unknown )
3651     {
3652         return false;
3653     }
3654 
3655     CPLAssert( dt.GetClass() != GEDTC_STRING );
3656     const bool bFastPath =
3657         ((m_dims.size() == 1 &&
3658           bufferStride[0] == 1 ) ||
3659          m_dims.empty() ) &&
3660         m_bPerfectDataTypeMatch &&
3661         bufferDataType == dt &&
3662         dt.GetSize() > 0;
3663     if( bFastPath )
3664     {
3665         int ret = nc_put_att(m_gid, m_varid, GetName().c_str(),
3666                              m_nAttType,
3667                              m_dims.empty() ? 1 : count[0],
3668                              pSrcBuffer);
3669         NCDF_ERR(ret);
3670         return ret == NC_NOERR;
3671     }
3672 
3673     const auto nElementSize = GetNCTypeSize(dt,
3674                                             m_bPerfectDataTypeMatch,
3675                                             m_nAttType);
3676     if( nElementSize == 0 )
3677         return false;
3678     const auto nInputDTSize = bufferDataType.GetSize();
3679     std::vector<GByte> abyBuffer(
3680         static_cast<size_t>(GetTotalElementsCount()) * nElementSize);
3681 
3682     const GByte* pabySrcBuffer = static_cast<const GByte*>(pSrcBuffer);
3683     auto pabyDstBuffer = &abyBuffer[0];
3684     for(size_t i = 0; i < (m_dims.empty() ? 1 : count[0]); i++ )
3685     {
3686         if( !m_bPerfectDataTypeMatch )
3687         {
3688             if( m_nAttType == NC_BYTE )
3689             {
3690                 short s;
3691                 GDALExtendedDataType::CopyValue(pabySrcBuffer, bufferDataType,
3692                             &s, dt);
3693                 signed char c = static_cast<signed char>(s);
3694                 memcpy(pabyDstBuffer, &c, sizeof(c));
3695             }
3696             else if( m_nAttType == NC_INT64 )
3697             {
3698                 double d;
3699                 GDALExtendedDataType::CopyValue(pabySrcBuffer, bufferDataType,
3700                             &d, dt);
3701                 GInt64 v = static_cast<GInt64>(d);
3702                 memcpy(pabyDstBuffer, &v, sizeof(v));
3703             }
3704             else if( m_nAttType == NC_UINT64 )
3705             {
3706                 double d;
3707                 GDALExtendedDataType::CopyValue(pabySrcBuffer, bufferDataType,
3708                             &d, dt);
3709                 GUInt64 v = static_cast<GUInt64>(d);
3710                 memcpy(pabyDstBuffer, &v, sizeof(v));
3711             }
3712             else
3713             {
3714                 CPLAssert(false);
3715             }
3716         }
3717         else
3718         {
3719             GDALExtendedDataType::CopyValue(pabySrcBuffer, bufferDataType,
3720                         pabyDstBuffer, dt);
3721         }
3722 
3723         if( !m_dims.empty() )
3724         {
3725             pabySrcBuffer += nInputDTSize * bufferStride[0];
3726             pabyDstBuffer += nElementSize;
3727         }
3728     }
3729 
3730     int ret = nc_put_att(m_gid, m_varid, GetName().c_str(),
3731                             m_nAttType,
3732                             m_dims.empty() ? 1 : count[0],
3733                             &abyBuffer[0]);
3734     NCDF_ERR(ret);
3735     return ret == NC_NOERR;
3736 }
3737 
3738 /************************************************************************/
3739 /*                           OpenMultiDim()                             */
3740 /************************************************************************/
3741 
OpenMultiDim(GDALOpenInfo * poOpenInfo)3742 GDALDataset *netCDFDataset::OpenMultiDim( GDALOpenInfo *poOpenInfo )
3743 {
3744 
3745     CPLMutexHolderD(&hNCMutex);
3746 
3747     CPLReleaseMutex(hNCMutex);  // Release mutex otherwise we'll deadlock with
3748                                 // GDALDataset own mutex.
3749     netCDFDataset *poDS = new netCDFDataset();
3750     CPLAcquireMutex(hNCMutex, 1000.0);
3751 
3752     auto poSharedResources(std::make_shared<netCDFSharedResources>());
3753 
3754     // For example to open DAP datasets
3755     if( STARTS_WITH_CI(poOpenInfo->pszFilename, "NETCDF:") )
3756     {
3757         poSharedResources->m_osFilename = poOpenInfo->pszFilename + strlen("NETCDF:");
3758         if( !poSharedResources->m_osFilename.empty() &&
3759             poSharedResources->m_osFilename[0] == '"' &&
3760             poSharedResources->m_osFilename.back() == '"' )
3761         {
3762             poSharedResources->m_osFilename = poSharedResources->m_osFilename.
3763                 substr(1, poSharedResources->m_osFilename.size() - 2);
3764         }
3765     }
3766     else
3767         poSharedResources->m_osFilename = poOpenInfo->pszFilename;
3768 
3769     poDS->SetDescription(poOpenInfo->pszFilename);
3770     poDS->papszOpenOptions = CSLDuplicate(poOpenInfo->papszOpenOptions);
3771 
3772 #ifdef ENABLE_NCDUMP
3773     const char* pszHeader =
3774                 reinterpret_cast<const char*>(poOpenInfo->pabyHeader);
3775     if( poOpenInfo->fpL != nullptr &&
3776         STARTS_WITH(pszHeader, "netcdf ") &&
3777         strstr(pszHeader, "dimensions:") &&
3778         strstr(pszHeader, "variables:") )
3779     {
3780         // By default create a temporary file that will be destroyed,
3781         // unless NETCDF_TMP_FILE is defined. Can be useful to see which
3782         // netCDF file has been generated from a potential fuzzed input.
3783         poSharedResources->m_osFilename = CPLGetConfigOption("NETCDF_TMP_FILE", "");
3784         if( poSharedResources->m_osFilename.empty() )
3785         {
3786             poSharedResources->m_bFileToDestroyAtClosing = true;
3787             poSharedResources->m_osFilename = CPLGenerateTempFilename("netcdf_tmp");
3788         }
3789         if( !netCDFDatasetCreateTempFile( NCDF_FORMAT_NC4,
3790                                           poSharedResources->m_osFilename,
3791                                           poOpenInfo->fpL ) )
3792         {
3793             CPLReleaseMutex(hNCMutex);  // Release mutex otherwise we'll deadlock
3794                                         // with GDALDataset own mutex.
3795             delete poDS;
3796             CPLAcquireMutex(hNCMutex, 1000.0);
3797             return nullptr;
3798         }
3799         poDS->eFormat = NCDF_FORMAT_NC4;
3800     }
3801 #endif
3802 
3803     // Try opening the dataset.
3804 #if defined(NCDF_DEBUG) && defined(ENABLE_UFFD)
3805     CPLDebug("GDAL_netCDF", "calling nc_open_mem(%s)", poSharedResources->m_osFilename.c_str());
3806 #elseif defined(NCDF_DEBUG) && !defined(ENABLE_UFFD)
3807     CPLDebug("GDAL_netCDF", "calling nc_open(%s)", poSharedResources->m_osFilename.c_str());
3808 #endif
3809     int cdfid = -1;
3810     const int nMode = (poOpenInfo->nOpenFlags & GDAL_OF_UPDATE) != 0 ? NC_WRITE : NC_NOWRITE;
3811     CPLString osFilenameForNCOpen(poSharedResources->m_osFilename);
3812 #ifdef WIN32
3813     if( CPLTestBool(CPLGetConfigOption( "GDAL_FILENAME_IS_UTF8", "YES" ) ) )
3814     {
3815         char* pszTemp = CPLRecode( osFilenameForNCOpen, CPL_ENC_UTF8, "CP_ACP" );
3816         osFilenameForNCOpen = pszTemp;
3817         CPLFree(pszTemp);
3818     }
3819 #endif
3820     int status2;
3821 
3822 #ifdef ENABLE_UFFD
3823     bool bVsiFile = !strncmp(osFilenameForNCOpen, "/vsi", strlen("/vsi"));
3824     bool bReadOnly = (poOpenInfo->eAccess == GA_ReadOnly);
3825     void * pVma = nullptr;
3826     uint64_t nVmaSize = 0;
3827     cpl_uffd_context * pCtx = nullptr;
3828 
3829     if ( bVsiFile && bReadOnly && CPLIsUserFaultMappingSupported() )
3830       pCtx = CPLCreateUserFaultMapping(osFilenameForNCOpen, &pVma, &nVmaSize);
3831     if (pCtx != nullptr && pVma != nullptr && nVmaSize > 0)
3832       status2 = nc_open_mem(osFilenameForNCOpen, nMode, static_cast<size_t>(nVmaSize), pVma, &cdfid);
3833     else
3834       status2 = nc_open(osFilenameForNCOpen, nMode, &cdfid);
3835     poSharedResources->m_pUffdCtx = pCtx;
3836 #else
3837     status2 = nc_open(osFilenameForNCOpen, nMode, &cdfid);
3838 #endif
3839     if( status2 != NC_NOERR )
3840     {
3841 #ifdef NCDF_DEBUG
3842         CPLDebug("GDAL_netCDF", "error opening");
3843 #endif
3844         CPLReleaseMutex(hNCMutex);  // Release mutex otherwise we'll deadlock
3845                                     // with GDALDataset own mutex.
3846         delete poDS;
3847         CPLAcquireMutex(hNCMutex, 1000.0);
3848         return nullptr;
3849     }
3850 #ifdef NCDF_DEBUG
3851     CPLDebug("GDAL_netCDF", "got cdfid=%d", cdfid);
3852 #endif
3853 
3854 #if defined(ENABLE_NCDUMP) && !defined(WIN32)
3855     // Try to destroy the temporary file right now on Unix
3856     if( poSharedResources->m_bFileToDestroyAtClosing )
3857     {
3858         if( VSIUnlink( poSharedResources->m_osFilename ) == 0 )
3859         {
3860             poSharedResources->m_bFileToDestroyAtClosing = false;
3861         }
3862     }
3863 #endif
3864     poSharedResources->m_bReadOnly = nMode == NC_NOWRITE;
3865     poSharedResources->m_cdfid = cdfid;
3866 
3867     // Is this a real netCDF file?
3868     int ndims;
3869     int ngatts;
3870     int nvars;
3871     int unlimdimid;
3872     int status = nc_inq(cdfid, &ndims, &nvars, &ngatts, &unlimdimid);
3873     if( status != NC_NOERR )
3874     {
3875         CPLReleaseMutex(hNCMutex);  // Release mutex otherwise we'll deadlock
3876                                     // with GDALDataset own mutex.
3877         delete poDS;
3878         CPLAcquireMutex(hNCMutex, 1000.0);
3879         return nullptr;
3880     }
3881 
3882     poDS->m_poRootGroup.reset(new netCDFGroup(poSharedResources, cdfid));
3883 
3884     poDS->TryLoadXML();
3885 
3886     return poDS;
3887 }
3888 
3889 /************************************************************************/
3890 /*                          GetRootGroup()                              */
3891 /************************************************************************/
3892 
GetRootGroup() const3893 std::shared_ptr<GDALGroup> netCDFDataset::GetRootGroup() const
3894 {
3895     return m_poRootGroup;
3896 }
3897 
3898 /************************************************************************/
3899 /*                      CreateMultiDimensional()                        */
3900 /************************************************************************/
3901 
CreateMultiDimensional(const char * pszFilename,CSLConstList,CSLConstList papszOptions)3902 GDALDataset* netCDFDataset::CreateMultiDimensional( const char * pszFilename,
3903                                                     CSLConstList /* papszRootGroupOptions */,
3904                                                     CSLConstList papszOptions )
3905 {
3906     CPLMutexHolderD(&hNCMutex);
3907 
3908     CPLReleaseMutex(hNCMutex);  // Release mutex otherwise we'll deadlock with
3909                                 // GDALDataset own mutex.
3910     netCDFDataset *poDS = new netCDFDataset();
3911     CPLAcquireMutex(hNCMutex, 1000.0);
3912     poDS->eAccess = GA_Update;
3913     poDS->osFilename = pszFilename;
3914 
3915     // process options.
3916     poDS->papszCreationOptions = CSLDuplicate(papszOptions);
3917     if( CSLFetchNameValue(papszOptions, "FORMAT") == nullptr )
3918     {
3919         poDS->papszCreationOptions = CSLSetNameValue(
3920             poDS->papszCreationOptions, "FORMAT", "NC4");
3921     }
3922     poDS->ProcessCreationOptions();
3923 
3924     // Create the dataset.
3925     CPLString osFilenameForNCCreate(pszFilename);
3926 #ifdef WIN32
3927     if( CPLTestBool(CPLGetConfigOption( "GDAL_FILENAME_IS_UTF8", "YES" ) ) )
3928     {
3929         char* pszTemp = CPLRecode( osFilenameForNCCreate, CPL_ENC_UTF8, "CP_ACP" );
3930         osFilenameForNCCreate = pszTemp;
3931         CPLFree(pszTemp);
3932     }
3933 #endif
3934     int cdfid = 0;
3935     int status = nc_create(osFilenameForNCCreate, poDS->nCreateMode, &cdfid);
3936     if( status != NC_NOERR )
3937     {
3938         CPLError(CE_Failure, CPLE_OpenFailed,
3939                  "Unable to create netCDF file %s (Error code %d): %s .",
3940                  pszFilename, status, nc_strerror(status));
3941         CPLReleaseMutex(hNCMutex);  // Release mutex otherwise we'll deadlock
3942                                     // with GDALDataset own mutex.
3943         delete poDS;
3944         CPLAcquireMutex(hNCMutex, 1000.0);
3945         return nullptr;
3946     }
3947 
3948     auto poSharedResources(std::make_shared<netCDFSharedResources>());
3949     poSharedResources->m_cdfid = cdfid;
3950     poSharedResources->m_bReadOnly = false;
3951     poSharedResources->m_bDefineMode = true;
3952     poDS->m_poRootGroup.reset(new netCDFGroup(poSharedResources, cdfid));
3953     const char* pszConventions = CSLFetchNameValueDef(
3954         papszOptions, "CONVENTIONS", NCDF_CONVENTIONS_CF_V1_6);
3955     if( !EQUAL(pszConventions, "") )
3956     {
3957         auto poAttr = poDS->m_poRootGroup->CreateAttribute(
3958             NCDF_CONVENTIONS, {}, GDALExtendedDataType::CreateString());
3959         if( poAttr )
3960             poAttr->Write(pszConventions);
3961     }
3962 
3963     return poDS;
3964 }
3965 
3966 #endif // NETCDF_HAS_NC4
3967