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