Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: netCDF read/write Driver
4 : * Purpose: Sentinel 3 SRAL/MWR Level 2 products
5 : * Author: Even Rouault <even.rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2019, Even Rouault <even.rouault at spatialys.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : // Example product:
14 : // https://scihub.copernicus.eu/dhus/odata/v1/Products('65b615b0-0db9-4ced-8020-eb17818f0c26')/$value
15 : // Specification:
16 : // https://sentinel.esa.int/documents/247904/2753172/Sentinel-3-Product-Data-Format-Specification-Level-2-Land
17 :
18 : #include "netcdfdataset.h"
19 :
20 : #include <limits>
21 :
22 : /************************************************************************/
23 : /* Sentinel3_SRAL_MWR_Layer */
24 : /************************************************************************/
25 :
26 : class Sentinel3_SRAL_MWR_Layer final : public OGRLayer
27 : {
28 : OGRFeatureDefn *m_poFDefn = nullptr;
29 : int m_cdfid;
30 : size_t m_nCurIdx = 0;
31 : size_t m_nFeatureCount = 0;
32 : CPLStringList m_aosMetadata{};
33 :
34 : struct VariableInfo
35 : {
36 : int varid;
37 : nc_type nctype;
38 : double scale;
39 : double offset;
40 : double nodata;
41 : };
42 :
43 : std::vector<VariableInfo> m_asVarInfo{};
44 : int m_iLongitude = -1;
45 : int m_iLatitude = -1;
46 : double m_dfLongScale = 1.0;
47 : double m_dfLongOffset = 0.0;
48 : double m_dfLatScale = 1.0;
49 : double m_dfLatOffset = 0.0;
50 :
51 : OGRFeature *TranslateFeature(size_t nIndex);
52 : OGRFeature *GetNextRawFeature();
53 :
54 : public:
55 : Sentinel3_SRAL_MWR_Layer(const std::string &name, int cdfid, int dimid);
56 : ~Sentinel3_SRAL_MWR_Layer();
57 :
58 10 : OGRFeatureDefn *GetLayerDefn() override
59 : {
60 10 : return m_poFDefn;
61 : }
62 :
63 : void ResetReading() override;
64 : OGRFeature *GetNextFeature() override;
65 : OGRFeature *GetFeature(GIntBig nFID) override;
66 : GIntBig GetFeatureCount(int bForce) override;
67 : int TestCapability(const char *pszCap) override;
68 : char **GetMetadata(const char *pszDomain) override;
69 : const char *GetMetadataItem(const char *pszKey,
70 : const char *pszDomain) override;
71 : };
72 :
73 : /************************************************************************/
74 : /* Sentinel3_SRAL_MWR_Layer() */
75 : /************************************************************************/
76 :
77 3 : Sentinel3_SRAL_MWR_Layer::Sentinel3_SRAL_MWR_Layer(const std::string &name,
78 3 : int cdfid, int dimid)
79 3 : : m_cdfid(cdfid)
80 : {
81 3 : m_poFDefn = new OGRFeatureDefn(name.c_str());
82 3 : m_poFDefn->SetGeomType(wkbPoint);
83 3 : m_poFDefn->Reference();
84 3 : SetDescription(name.c_str());
85 :
86 3 : nc_inq_dimlen(cdfid, dimid, &m_nFeatureCount);
87 :
88 3 : int nVars = 0;
89 3 : NCDF_ERR(nc_inq(cdfid, nullptr, &nVars, nullptr, nullptr));
90 45 : for (int iVar = 0; iVar < nVars; iVar++)
91 : {
92 42 : int nVarDims = 0;
93 42 : NCDF_ERR(nc_inq_varndims(cdfid, iVar, &nVarDims));
94 42 : if (nVarDims != 1)
95 34 : continue;
96 :
97 42 : int vardimid = -1;
98 42 : NCDF_ERR(nc_inq_vardimid(cdfid, iVar, &vardimid));
99 42 : if (vardimid != dimid)
100 28 : continue;
101 :
102 14 : char szVarName[NC_MAX_NAME + 1] = {};
103 14 : NCDF_ERR(nc_inq_varname(cdfid, iVar, szVarName));
104 :
105 14 : nc_type vartype = NC_NAT;
106 14 : nc_inq_vartype(cdfid, iVar, &vartype);
107 :
108 14 : int nbAttr = 0;
109 14 : NCDF_ERR(nc_inq_varnatts(cdfid, iVar, &nbAttr));
110 14 : std::string scaleFactor;
111 14 : std::string offset;
112 14 : std::string fillValue;
113 14 : CPLStringList aosMetadata;
114 93 : for (int iAttr = 0; iAttr < nbAttr; iAttr++)
115 : {
116 : char szAttrName[NC_MAX_NAME + 1];
117 79 : szAttrName[0] = 0;
118 79 : NCDF_ERR(nc_inq_attname(cdfid, iVar, iAttr, szAttrName));
119 79 : char *pszMetaTemp = nullptr;
120 158 : if (NCDFGetAttr(cdfid, iVar, szAttrName, &pszMetaTemp) == CE_None &&
121 79 : pszMetaTemp)
122 : {
123 79 : if (EQUAL(szAttrName, "scale_factor"))
124 : {
125 9 : scaleFactor = pszMetaTemp;
126 : }
127 70 : else if (EQUAL(szAttrName, "add_offset"))
128 : {
129 9 : offset = pszMetaTemp;
130 : }
131 61 : else if (EQUAL(szAttrName, "_FillValue"))
132 : {
133 5 : fillValue = pszMetaTemp;
134 : }
135 56 : else if (!EQUAL(szAttrName, "coordinates"))
136 : {
137 51 : aosMetadata.SetNameValue(szAttrName, pszMetaTemp);
138 : }
139 : }
140 79 : CPLFree(pszMetaTemp);
141 : }
142 :
143 : const char *pszStandardName =
144 14 : aosMetadata.FetchNameValue("standard_name");
145 14 : if (pszStandardName)
146 : {
147 10 : if (EQUAL(pszStandardName, "longitude") && vartype == NC_INT)
148 : {
149 3 : m_iLongitude = iVar;
150 3 : if (!scaleFactor.empty())
151 3 : m_dfLongScale = CPLAtof(scaleFactor.c_str());
152 3 : if (!offset.empty())
153 3 : m_dfLongOffset = CPLAtof(offset.c_str());
154 3 : continue;
155 : }
156 7 : if (EQUAL(pszStandardName, "latitude") && vartype == NC_INT)
157 : {
158 3 : m_iLatitude = iVar;
159 3 : if (!scaleFactor.empty())
160 3 : m_dfLatScale = CPLAtof(scaleFactor.c_str());
161 3 : if (!offset.empty())
162 3 : m_dfLatOffset = CPLAtof(offset.c_str());
163 3 : continue;
164 : }
165 : }
166 :
167 35 : for (int i = 0; i < aosMetadata.size(); i++)
168 : m_aosMetadata.AddString(
169 27 : (std::string(szVarName) + '_' + aosMetadata[i]).c_str());
170 :
171 8 : OGRFieldType eType = OFTReal;
172 8 : if (!scaleFactor.empty())
173 : {
174 : // do nothing
175 : }
176 5 : else if (!offset.empty())
177 : {
178 : // do nothing
179 : }
180 5 : else if (vartype == NC_BYTE || vartype == NC_SHORT ||
181 4 : vartype == NC_INT || vartype == NC_USHORT ||
182 3 : vartype == NC_UINT)
183 : {
184 2 : eType = OFTInteger;
185 : }
186 16 : OGRFieldDefn oField(szVarName, eType);
187 8 : m_poFDefn->AddFieldDefn(&oField);
188 : VariableInfo varInfo;
189 8 : varInfo.varid = iVar;
190 8 : varInfo.nctype = vartype;
191 8 : varInfo.scale =
192 8 : scaleFactor.empty() ? 1.0 : CPLAtof(scaleFactor.c_str());
193 8 : varInfo.offset = offset.empty() ? 0.0 : CPLAtof(offset.c_str());
194 8 : varInfo.nodata = fillValue.empty()
195 8 : ? std::numeric_limits<double>::quiet_NaN()
196 5 : : CPLAtof(fillValue.c_str());
197 8 : m_asVarInfo.emplace_back(varInfo);
198 : }
199 3 : }
200 :
201 : /************************************************************************/
202 : /* ~Sentinel3_SRAL_MWR_Layer() */
203 : /************************************************************************/
204 :
205 6 : Sentinel3_SRAL_MWR_Layer::~Sentinel3_SRAL_MWR_Layer()
206 : {
207 3 : m_poFDefn->Release();
208 6 : }
209 :
210 : /************************************************************************/
211 : /* GetMetadata() */
212 : /************************************************************************/
213 :
214 1 : char **Sentinel3_SRAL_MWR_Layer::GetMetadata(const char *pszDomain)
215 : {
216 1 : if (pszDomain == nullptr || EQUAL(pszDomain, ""))
217 1 : return m_aosMetadata.List();
218 0 : return OGRLayer::GetMetadata(pszDomain);
219 : }
220 :
221 : /************************************************************************/
222 : /* GetMetadataItem() */
223 : /************************************************************************/
224 :
225 1 : const char *Sentinel3_SRAL_MWR_Layer::GetMetadataItem(const char *pszKey,
226 : const char *pszDomain)
227 : {
228 1 : if (pszDomain == nullptr || EQUAL(pszDomain, ""))
229 1 : return m_aosMetadata.FetchNameValue(pszKey);
230 0 : return OGRLayer::GetMetadataItem(pszKey, pszDomain);
231 : }
232 :
233 : /************************************************************************/
234 : /* ResetReading() */
235 : /************************************************************************/
236 :
237 8 : void Sentinel3_SRAL_MWR_Layer::ResetReading()
238 : {
239 8 : m_nCurIdx = 0;
240 8 : }
241 :
242 : /************************************************************************/
243 : /* GetFeatureCount() */
244 : /************************************************************************/
245 :
246 2 : GIntBig Sentinel3_SRAL_MWR_Layer::GetFeatureCount(int bForce)
247 : {
248 2 : if (m_poFilterGeom == nullptr && m_poAttrQuery == nullptr)
249 1 : return m_nFeatureCount;
250 1 : return OGRLayer::GetFeatureCount(bForce);
251 : }
252 :
253 : /************************************************************************/
254 : /* TestCapability() */
255 : /************************************************************************/
256 :
257 4 : int Sentinel3_SRAL_MWR_Layer::TestCapability(const char *pszCap)
258 : {
259 4 : if (EQUAL(pszCap, OLCFastFeatureCount))
260 1 : return m_poFilterGeom == nullptr && m_poAttrQuery == nullptr;
261 3 : if (EQUAL(pszCap, OLCRandomRead))
262 1 : return true;
263 2 : return false;
264 : }
265 :
266 : /************************************************************************/
267 : /* TranslateFeature() */
268 : /************************************************************************/
269 :
270 12 : OGRFeature *Sentinel3_SRAL_MWR_Layer::TranslateFeature(size_t nIndex)
271 : {
272 12 : OGRFeature *poFeat = new OGRFeature(m_poFDefn);
273 12 : poFeat->SetFID(nIndex + 1);
274 12 : if (m_iLongitude >= 0 && m_iLatitude >= 0)
275 : {
276 12 : int nLong = 0;
277 12 : int status = nc_get_var1_int(m_cdfid, m_iLongitude, &nIndex, &nLong);
278 12 : if (status == NC_NOERR)
279 : {
280 12 : int nLat = 0;
281 12 : status = nc_get_var1_int(m_cdfid, m_iLatitude, &nIndex, &nLat);
282 12 : if (status == NC_NOERR)
283 : {
284 12 : const double dfLong = nLong * m_dfLongScale + m_dfLongOffset;
285 12 : const double dfLat = nLat * m_dfLatScale + m_dfLatOffset;
286 12 : auto poGeom = new OGRPoint(dfLong, dfLat);
287 12 : auto poGeomField = m_poFDefn->GetGeomFieldDefn(0);
288 12 : poGeom->assignSpatialReference(poGeomField->GetSpatialRef());
289 12 : poFeat->SetGeometryDirectly(poGeom);
290 : }
291 : }
292 : }
293 :
294 66 : for (int i = 0; i < static_cast<int>(m_asVarInfo.size()); i++)
295 : {
296 54 : if (m_asVarInfo[i].nctype == NC_BYTE)
297 : {
298 10 : signed char nVal = 0;
299 10 : int status = nc_get_var1_schar(m_cdfid, m_asVarInfo[i].varid,
300 : &nIndex, &nVal);
301 10 : if (status == NC_NOERR && nVal != m_asVarInfo[i].nodata)
302 : {
303 4 : poFeat->SetField(i, nVal * m_asVarInfo[i].scale +
304 4 : m_asVarInfo[i].offset);
305 : }
306 : }
307 44 : else if (m_asVarInfo[i].nctype == NC_SHORT)
308 : {
309 10 : short nVal = 0;
310 10 : int status = nc_get_var1_short(m_cdfid, m_asVarInfo[i].varid,
311 : &nIndex, &nVal);
312 10 : if (status == NC_NOERR && nVal != m_asVarInfo[i].nodata)
313 : {
314 4 : poFeat->SetField(i, nVal * m_asVarInfo[i].scale +
315 4 : m_asVarInfo[i].offset);
316 : }
317 : }
318 34 : else if (m_asVarInfo[i].nctype == NC_USHORT)
319 : {
320 2 : unsigned short nVal = 0;
321 2 : int status = nc_get_var1_ushort(m_cdfid, m_asVarInfo[i].varid,
322 : &nIndex, &nVal);
323 2 : if (status == NC_NOERR && nVal != m_asVarInfo[i].nodata)
324 : {
325 1 : poFeat->SetField(i, nVal * m_asVarInfo[i].scale +
326 1 : m_asVarInfo[i].offset);
327 : }
328 : }
329 32 : else if (m_asVarInfo[i].nctype == NC_INT)
330 : {
331 10 : int nVal = 0;
332 : int status =
333 10 : nc_get_var1_int(m_cdfid, m_asVarInfo[i].varid, &nIndex, &nVal);
334 10 : if (status == NC_NOERR && nVal != m_asVarInfo[i].nodata)
335 : {
336 4 : poFeat->SetField(i, nVal * m_asVarInfo[i].scale +
337 4 : m_asVarInfo[i].offset);
338 : }
339 : }
340 22 : else if (m_asVarInfo[i].nctype == NC_UINT)
341 : {
342 10 : unsigned int nVal = 0;
343 : int status =
344 10 : nc_get_var1_uint(m_cdfid, m_asVarInfo[i].varid, &nIndex, &nVal);
345 10 : if (status == NC_NOERR && nVal != m_asVarInfo[i].nodata)
346 : {
347 4 : poFeat->SetField(i, nVal * m_asVarInfo[i].scale +
348 4 : m_asVarInfo[i].offset);
349 : }
350 : }
351 12 : else if (m_asVarInfo[i].nctype == NC_DOUBLE)
352 : {
353 12 : double dfVal = 0.0;
354 12 : int status = nc_get_var1_double(m_cdfid, m_asVarInfo[i].varid,
355 : &nIndex, &dfVal);
356 12 : if (status == NC_NOERR && dfVal != m_asVarInfo[i].nodata)
357 : {
358 12 : poFeat->SetField(i, dfVal * m_asVarInfo[i].scale +
359 12 : m_asVarInfo[i].offset);
360 : }
361 : }
362 : else
363 : {
364 0 : CPLDebug("netCDF", "Unhandled data type %d for %s",
365 0 : m_asVarInfo[i].nctype,
366 0 : m_poFDefn->GetFieldDefn(i)->GetNameRef());
367 : }
368 : }
369 :
370 12 : return poFeat;
371 : }
372 :
373 : /************************************************************************/
374 : /* GetNextRawFeature() */
375 : /************************************************************************/
376 :
377 16 : OGRFeature *Sentinel3_SRAL_MWR_Layer::GetNextRawFeature()
378 : {
379 16 : if (m_nCurIdx == m_nFeatureCount)
380 5 : return nullptr;
381 11 : OGRFeature *poFeat = TranslateFeature(m_nCurIdx);
382 11 : m_nCurIdx++;
383 11 : return poFeat;
384 : }
385 :
386 : /************************************************************************/
387 : /* GetFeature() */
388 : /************************************************************************/
389 :
390 3 : OGRFeature *Sentinel3_SRAL_MWR_Layer::GetFeature(GIntBig nFID)
391 : {
392 3 : if (nFID <= 0 || static_cast<size_t>(nFID) > m_nFeatureCount)
393 2 : return nullptr;
394 1 : return TranslateFeature(static_cast<size_t>(nFID - 1));
395 : }
396 :
397 : /************************************************************************/
398 : /* GetNextFeature() */
399 : /************************************************************************/
400 :
401 16 : OGRFeature *Sentinel3_SRAL_MWR_Layer::GetNextFeature()
402 : {
403 : while (true)
404 : {
405 16 : OGRFeature *poFeature = GetNextRawFeature();
406 16 : if (poFeature == nullptr)
407 5 : return nullptr;
408 :
409 26 : if ((m_poFilterGeom == nullptr ||
410 18 : FilterGeometry(poFeature->GetGeomFieldRef(m_iGeomFieldFilter))) &&
411 7 : (m_poAttrQuery == nullptr || m_poAttrQuery->Evaluate(poFeature)))
412 5 : return poFeature;
413 :
414 6 : delete poFeature;
415 6 : }
416 : }
417 :
418 : /************************************************************************/
419 : /* ProcessSentinel3_SRAL_MWR() */
420 : /************************************************************************/
421 :
422 1 : void netCDFDataset::ProcessSentinel3_SRAL_MWR()
423 : {
424 1 : int nDimCount = -1;
425 1 : int status = nc_inq_ndims(cdfid, &nDimCount);
426 1 : NCDF_ERR(status);
427 1 : if (status != NC_NOERR || nDimCount == 0 || nDimCount > 1000)
428 0 : return;
429 1 : std::vector<int> dimIds(nDimCount);
430 1 : int nDimCount2 = -1;
431 1 : status = nc_inq_dimids(cdfid, &nDimCount2, &dimIds[0], FALSE);
432 1 : NCDF_ERR(status);
433 1 : if (status != NC_NOERR)
434 0 : return;
435 1 : CPLAssert(nDimCount == nDimCount2);
436 :
437 1 : OGRSpatialReference *poSRS = nullptr;
438 : const char *pszSemiMajor =
439 1 : CSLFetchNameValue(papszMetadata, "NC_GLOBAL#semi_major_ellipsoid_axis");
440 : const char *pszFlattening =
441 1 : CSLFetchNameValue(papszMetadata, "NC_GLOBAL#ellipsoid_flattening");
442 2 : if (pszSemiMajor && EQUAL(pszSemiMajor, "6378137") && pszFlattening &&
443 1 : fabs(CPLAtof(pszFlattening) - 0.00335281066474748) < 1e-16)
444 : {
445 : int iAttr =
446 1 : CSLFindName(papszMetadata, "NC_GLOBAL#semi_major_ellipsoid_axis");
447 1 : if (iAttr >= 0)
448 1 : papszMetadata = CSLRemoveStrings(papszMetadata, iAttr, 1, nullptr);
449 1 : iAttr = CSLFindName(papszMetadata, "NC_GLOBAL#ellipsoid_flattening");
450 1 : if (iAttr >= 0)
451 1 : papszMetadata = CSLRemoveStrings(papszMetadata, iAttr, 1, nullptr);
452 1 : poSRS = new OGRSpatialReference();
453 1 : poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
454 1 : poSRS->importFromEPSG(4326);
455 : }
456 :
457 4 : for (int i = 0; i < nDimCount; i++)
458 : {
459 3 : char szDimName[NC_MAX_NAME + 1] = {};
460 3 : status = nc_inq_dimname(cdfid, dimIds[i], szDimName);
461 3 : NCDF_ERR(status);
462 3 : if (status != NC_NOERR)
463 0 : break;
464 6 : std::string name(CPLGetBasename(GetDescription()));
465 3 : name += '_';
466 3 : name += szDimName;
467 : std::shared_ptr<OGRLayer> poLayer(
468 9 : new Sentinel3_SRAL_MWR_Layer(name.c_str(), cdfid, dimIds[i]));
469 3 : auto poGeomField = poLayer->GetLayerDefn()->GetGeomFieldDefn(0);
470 3 : if (poGeomField)
471 3 : poGeomField->SetSpatialRef(poSRS);
472 3 : papoLayers.emplace_back(poLayer);
473 : }
474 :
475 1 : if (poSRS)
476 1 : poSRS->Release();
477 : }
|
