Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GRIB Driver
4 : * Purpose: GDALDataset driver for GRIB translator for read support
5 : * Author: Bas Retsios, retsios@itc.nl
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2007, ITC
9 : * Copyright (c) 2008-2017, Even Rouault <even dot rouault at spatialys dot com>
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ******************************************************************************
29 : *
30 : */
31 :
32 : #include "cpl_port.h"
33 : #include "gribdataset.h"
34 : #include "gribdrivercore.h"
35 :
36 : #include <cerrno>
37 : #include <cmath>
38 : #include <cstddef>
39 : #include <cstdio>
40 : #include <cstdlib>
41 : #include <cstring>
42 : #if HAVE_FCNTL_H
43 : #include <fcntl.h>
44 : #endif
45 :
46 : #include <algorithm>
47 : #include <set>
48 : #include <string>
49 : #include <vector>
50 :
51 : #include "cpl_conv.h"
52 : #include "cpl_error.h"
53 : #include "cpl_multiproc.h"
54 : #include "cpl_string.h"
55 : #include "cpl_vsi.h"
56 : #include "cpl_time.h"
57 : #include "degrib/degrib/degrib2.h"
58 : #include "degrib/degrib/inventory.h"
59 : #include "degrib/degrib/meta.h"
60 : #include "degrib/degrib/metaname.h"
61 : #include "degrib/degrib/myerror.h"
62 : #include "degrib/degrib/type.h"
63 : CPL_C_START
64 : #include "degrib/g2clib/grib2.h"
65 : #include "degrib/g2clib/pdstemplates.h"
66 : CPL_C_END
67 : #include "gdal.h"
68 : #include "gdal_frmts.h"
69 : #include "gdal_pam.h"
70 : #include "gdal_priv.h"
71 : #include "ogr_spatialref.h"
72 : #include "memdataset.h"
73 :
74 : static CPLMutex *hGRIBMutex = nullptr;
75 :
76 : /************************************************************************/
77 : /* ConvertUnitInText() */
78 : /************************************************************************/
79 :
80 908 : static CPLString ConvertUnitInText(bool bMetricUnits, const char *pszTxt)
81 : {
82 908 : if (pszTxt == nullptr)
83 0 : return CPLString();
84 908 : if (!bMetricUnits)
85 4 : return pszTxt;
86 :
87 1808 : CPLString osRes(pszTxt);
88 904 : size_t iPos = osRes.find("[K]");
89 904 : if (iPos != std::string::npos)
90 96 : osRes = osRes.substr(0, iPos) + "[C]" + osRes.substr(iPos + 3);
91 904 : return osRes;
92 : }
93 :
94 : /************************************************************************/
95 : /* Lon360to180() */
96 : /************************************************************************/
97 :
98 355 : static inline double Lon360to180(double lon)
99 : {
100 355 : if (lon == 180)
101 0 : return 180;
102 355 : return fmod(lon + 180, 360) - 180;
103 : }
104 :
105 : /************************************************************************/
106 : /* GRIBRasterBand() */
107 : /************************************************************************/
108 :
109 431 : GRIBRasterBand::GRIBRasterBand(GRIBDataset *poDSIn, int nBandIn,
110 431 : inventoryType *psInv)
111 431 : : start(psInv->start), subgNum(psInv->subgNum),
112 862 : longFstLevel(CPLStrdup(psInv->longFstLevel)), m_Grib_Data(nullptr),
113 431 : m_Grib_MetaData(nullptr), nGribDataXSize(poDSIn->nRasterXSize),
114 431 : nGribDataYSize(poDSIn->nRasterYSize), m_nGribVersion(psInv->GribVersion),
115 431 : m_bHasLookedForNoData(false), m_dfNoData(0.0), m_bHasNoData(false)
116 :
117 : {
118 431 : poDS = poDSIn;
119 431 : nBand = nBandIn;
120 :
121 : // Let user do -ot Float32 if needed for saving space, GRIB contains
122 : // Float64 (though not fully utilized most of the time).
123 431 : eDataType = GDT_Float64;
124 :
125 431 : nBlockXSize = poDSIn->nRasterXSize;
126 431 : nBlockYSize = 1;
127 :
128 431 : if (psInv->unitName != nullptr && psInv->comment != nullptr &&
129 405 : psInv->element != nullptr)
130 : {
131 405 : bLoadedMetadata = true;
132 : const char *pszGribNormalizeUnits =
133 405 : CPLGetConfigOption("GRIB_NORMALIZE_UNITS", "YES");
134 405 : bool bMetricUnits = CPLTestBool(pszGribNormalizeUnits);
135 :
136 405 : SetMetadataItem("GRIB_UNIT",
137 810 : ConvertUnitInText(bMetricUnits, psInv->unitName));
138 405 : SetMetadataItem("GRIB_COMMENT",
139 810 : ConvertUnitInText(bMetricUnits, psInv->comment));
140 405 : SetMetadataItem("GRIB_ELEMENT", psInv->element);
141 405 : SetMetadataItem("GRIB_SHORT_NAME", psInv->shortFstLevel);
142 405 : SetMetadataItem("GRIB_REF_TIME",
143 810 : CPLString().Printf("%.0f", psInv->refTime));
144 405 : SetMetadataItem("GRIB_VALID_TIME",
145 810 : CPLString().Printf("%.0f", psInv->validTime));
146 405 : SetMetadataItem("GRIB_FORECAST_SECONDS",
147 810 : CPLString().Printf("%.0f", psInv->foreSec));
148 : }
149 431 : }
150 :
151 : /************************************************************************/
152 : /* FindMetaData() */
153 : /************************************************************************/
154 :
155 779 : void GRIBRasterBand::FindMetaData()
156 : {
157 779 : if (bLoadedMetadata)
158 753 : return;
159 26 : if (m_Grib_MetaData == nullptr)
160 : {
161 : grib_MetaData *metaData;
162 5 : GRIBDataset *poGDS = static_cast<GRIBDataset *>(poDS);
163 5 : GRIBRasterBand::ReadGribData(poGDS->fp, start, subgNum, nullptr,
164 : &metaData);
165 5 : if (metaData == nullptr)
166 0 : return;
167 5 : m_Grib_MetaData = metaData;
168 : }
169 26 : bLoadedMetadata = true;
170 26 : m_nGribVersion = m_Grib_MetaData->GribVersion;
171 :
172 : const char *pszGribNormalizeUnits =
173 26 : CPLGetConfigOption("GRIB_NORMALIZE_UNITS", "YES");
174 26 : bool bMetricUnits = CPLTestBool(pszGribNormalizeUnits);
175 :
176 26 : GDALRasterBand::SetMetadataItem(
177 : "GRIB_UNIT",
178 52 : ConvertUnitInText(bMetricUnits, m_Grib_MetaData->unitName));
179 26 : GDALRasterBand::SetMetadataItem(
180 : "GRIB_COMMENT",
181 52 : ConvertUnitInText(bMetricUnits, m_Grib_MetaData->comment));
182 :
183 26 : GDALRasterBand::SetMetadataItem("GRIB_ELEMENT", m_Grib_MetaData->element);
184 26 : GDALRasterBand::SetMetadataItem("GRIB_SHORT_NAME",
185 26 : m_Grib_MetaData->shortFstLevel);
186 :
187 26 : if (m_nGribVersion == 2)
188 : {
189 14 : GDALRasterBand::SetMetadataItem(
190 : "GRIB_REF_TIME",
191 28 : CPLString().Printf("%.0f", m_Grib_MetaData->pds2.refTime));
192 14 : GDALRasterBand::SetMetadataItem(
193 : "GRIB_VALID_TIME",
194 28 : CPLString().Printf("%.0f", m_Grib_MetaData->pds2.sect4.validTime));
195 : }
196 12 : else if (m_nGribVersion == 1)
197 : {
198 12 : GDALRasterBand::SetMetadataItem(
199 : "GRIB_REF_TIME",
200 24 : CPLString().Printf("%.0f", m_Grib_MetaData->pds1.refTime));
201 12 : GDALRasterBand::SetMetadataItem(
202 : "GRIB_VALID_TIME",
203 24 : CPLString().Printf("%.0f", m_Grib_MetaData->pds1.validTime));
204 : }
205 :
206 26 : GDALRasterBand::SetMetadataItem(
207 : "GRIB_FORECAST_SECONDS",
208 52 : CPLString().Printf("%d", m_Grib_MetaData->deltTime));
209 : }
210 :
211 : /************************************************************************/
212 : /* FindTrueStart() */
213 : /* */
214 : /* Scan after the official start of the message to find its */
215 : /* true starting offset. */
216 : /************************************************************************/
217 842 : vsi_l_offset GRIBRasterBand::FindTrueStart(VSILFILE *fp, vsi_l_offset start)
218 : {
219 : // GRIB messages can be preceded by "garbage". GRIB2Inventory()
220 : // does not return the offset to the real start of the message
221 : char szHeader[1024 + 1];
222 842 : VSIFSeekL(fp, start, SEEK_SET);
223 : const int nRead =
224 842 : static_cast<int>(VSIFReadL(szHeader, 1, sizeof(szHeader) - 1, fp));
225 842 : szHeader[nRead] = 0;
226 : // Find the real offset of the fist message
227 842 : int nOffsetFirstMessage = 0;
228 2524 : for (int j = 0; j + 3 < nRead; j++)
229 : {
230 2524 : if (STARTS_WITH_CI(szHeader + j, "GRIB")
231 : #ifdef ENABLE_TDLP
232 : || STARTS_WITH_CI(szHeader + j, "TDLP")
233 : #endif
234 : )
235 : {
236 842 : nOffsetFirstMessage = j;
237 842 : break;
238 : }
239 : }
240 842 : return start + nOffsetFirstMessage;
241 : }
242 :
243 : /************************************************************************/
244 : /* FindPDSTemplateGRIB2() */
245 : /* */
246 : /* Scan the file for the PDS template info and represent it as */
247 : /* metadata. */
248 : /************************************************************************/
249 :
250 883 : void GRIBRasterBand::FindPDSTemplateGRIB2()
251 :
252 : {
253 883 : CPLAssert(m_nGribVersion == 2);
254 :
255 883 : if (bLoadedPDS)
256 564 : return;
257 319 : bLoadedPDS = true;
258 :
259 319 : GRIBDataset *poGDS = static_cast<GRIBDataset *>(poDS);
260 319 : start = FindTrueStart(poGDS->fp, start);
261 :
262 : // Read section 0
263 : GByte abySection0[16];
264 638 : if (VSIFSeekL(poGDS->fp, start, SEEK_SET) != 0 ||
265 319 : VSIFReadL(abySection0, 16, 1, poGDS->fp) != 1)
266 : {
267 0 : CPLDebug("GRIB", "Cannot read leading bytes of section 0");
268 0 : return;
269 : }
270 319 : GByte nDiscipline = abySection0[7 - 1];
271 319 : CPLString osDiscipline;
272 319 : osDiscipline = CPLString().Printf("%d", nDiscipline);
273 : static const char *const table00[] = {"Meteorological",
274 : "Hydrological",
275 : "Land Surface",
276 : "Space products",
277 : "Space products",
278 : "Reserved",
279 : "Reserved",
280 : "Reserved",
281 : "Reserved",
282 : "Reserved",
283 : "Oceanographic Products"};
284 319 : m_nDisciplineCode = nDiscipline;
285 319 : if (nDiscipline < CPL_ARRAYSIZE(table00))
286 : {
287 318 : m_osDisciplineName = table00[nDiscipline];
288 636 : osDiscipline += CPLString("(") +
289 1272 : CPLString(table00[nDiscipline]).replaceAll(' ', '_') +
290 318 : ")";
291 : }
292 :
293 319 : GDALRasterBand::SetMetadataItem("GRIB_DISCIPLINE", osDiscipline.c_str());
294 :
295 319 : GByte abyHead[5] = {0};
296 :
297 319 : if (VSIFReadL(abyHead, 5, 1, poGDS->fp) != 1)
298 : {
299 0 : CPLDebug("GRIB", "Cannot read 5 leading bytes past section 0");
300 0 : return;
301 : }
302 :
303 319 : GUInt32 nSectSize = 0;
304 319 : if (abyHead[4] == 1)
305 : {
306 319 : memcpy(&nSectSize, abyHead, 4);
307 319 : CPL_MSBPTR32(&nSectSize);
308 319 : if (nSectSize >= 21 && nSectSize <= 100000 /* arbitrary upper limit */)
309 : {
310 319 : GByte *pabyBody = static_cast<GByte *>(CPLMalloc(nSectSize));
311 319 : memcpy(pabyBody, abyHead, 5);
312 319 : VSIFReadL(pabyBody + 5, 1, nSectSize - 5, poGDS->fp);
313 :
314 638 : CPLString osIDS;
315 319 : unsigned short nCenter = static_cast<unsigned short>(
316 319 : pabyBody[6 - 1] * 256 + pabyBody[7 - 1]);
317 319 : if (nCenter != GRIB2MISSING_u1 && nCenter != GRIB2MISSING_u2)
318 : {
319 86 : osIDS += "CENTER=";
320 86 : m_nCenter = nCenter;
321 86 : osIDS += CPLSPrintf("%d", nCenter);
322 86 : const char *pszCenter = centerLookup(nCenter);
323 86 : if (pszCenter)
324 : {
325 86 : m_osCenterName = pszCenter;
326 86 : osIDS += CPLString("(") + pszCenter + ")";
327 : }
328 : }
329 :
330 319 : unsigned short nSubCenter = static_cast<unsigned short>(
331 319 : pabyBody[8 - 1] * 256 + pabyBody[9 - 1]);
332 319 : if (nSubCenter != GRIB2MISSING_u2)
333 : {
334 112 : if (!osIDS.empty())
335 71 : osIDS += " ";
336 112 : osIDS += "SUBCENTER=";
337 112 : osIDS += CPLSPrintf("%d", nSubCenter);
338 112 : m_nSubCenter = nSubCenter;
339 112 : const char *pszSubCenter = subCenterLookup(nCenter, nSubCenter);
340 112 : if (pszSubCenter)
341 : {
342 3 : m_osSubCenterName = pszSubCenter;
343 3 : osIDS += CPLString("(") + pszSubCenter + ")";
344 : }
345 : }
346 :
347 319 : if (!osIDS.empty())
348 127 : osIDS += " ";
349 319 : osIDS += "MASTER_TABLE=";
350 319 : osIDS += CPLSPrintf("%d", pabyBody[10 - 1]);
351 319 : osIDS += " ";
352 319 : osIDS += "LOCAL_TABLE=";
353 319 : osIDS += CPLSPrintf("%d", pabyBody[11 - 1]);
354 319 : osIDS += " ";
355 319 : osIDS += "SIGNF_REF_TIME=";
356 319 : unsigned nSignRefTime = pabyBody[12 - 1];
357 319 : osIDS += CPLSPrintf("%d", nSignRefTime);
358 : static const char *const table12[] = {
359 : "Analysis", "Start of Forecast", "Verifying time of forecast",
360 : "Observation time"};
361 319 : if (nSignRefTime < CPL_ARRAYSIZE(table12))
362 : {
363 319 : m_osSignRefTimeName = table12[nSignRefTime];
364 638 : osIDS += CPLString("(") +
365 1276 : CPLString(table12[nSignRefTime]).replaceAll(' ', '_') +
366 319 : ")";
367 : }
368 319 : osIDS += " ";
369 319 : osIDS += "REF_TIME=";
370 : m_osRefTime =
371 : CPLSPrintf("%04d-%02d-%02dT%02d:%02d:%02dZ",
372 319 : pabyBody[13 - 1] * 256 + pabyBody[14 - 1],
373 319 : pabyBody[15 - 1], pabyBody[16 - 1], pabyBody[17 - 1],
374 319 : pabyBody[18 - 1], pabyBody[19 - 1]);
375 319 : osIDS += m_osRefTime;
376 319 : osIDS += " ";
377 319 : osIDS += "PROD_STATUS=";
378 319 : unsigned nProdStatus = pabyBody[20 - 1];
379 319 : osIDS += CPLSPrintf("%d", nProdStatus);
380 : static const char *const table13[] = {
381 : "Operational", "Operational test",
382 : "Research", "Re-analysis",
383 : "TIGGE", "TIGGE test",
384 : "S2S operational", "S2S test",
385 : "UERRA", "UERRA test"};
386 319 : if (nProdStatus < CPL_ARRAYSIZE(table13))
387 : {
388 126 : m_osProductionStatus = table13[nProdStatus];
389 252 : osIDS += CPLString("(") +
390 504 : CPLString(table13[nProdStatus]).replaceAll(' ', '_') +
391 126 : ")";
392 : }
393 319 : osIDS += " ";
394 319 : osIDS += "TYPE=";
395 319 : unsigned nType = pabyBody[21 - 1];
396 319 : osIDS += CPLSPrintf("%d", nType);
397 : static const char *const table14[] = {
398 : "Analysis",
399 : "Forecast",
400 : "Analysis and forecast",
401 : "Control forecast",
402 : "Perturbed forecast",
403 : "Control and perturbed forecast",
404 : "Processed satellite observations",
405 : "Processed radar observations",
406 : "Event Probability"};
407 319 : if (nType < CPL_ARRAYSIZE(table14))
408 : {
409 126 : m_osType = table14[nType];
410 252 : osIDS += CPLString("(") +
411 378 : CPLString(table14[nType]).replaceAll(' ', '_') + ")";
412 : }
413 :
414 319 : GDALRasterBand::SetMetadataItem("GRIB_IDS", osIDS);
415 :
416 319 : CPLFree(pabyBody);
417 : }
418 :
419 319 : if (VSIFReadL(abyHead, 5, 1, poGDS->fp) != 1)
420 : {
421 0 : CPLDebug("GRIB", "Cannot read 5 leading bytes past section 1");
422 0 : return;
423 : }
424 : }
425 :
426 319 : if (subgNum > 0)
427 : {
428 : // If we are a subgrid, then iterate over all preceding subgrids
429 14 : for (int iSubMessage = 0; iSubMessage < subgNum;)
430 : {
431 12 : memcpy(&nSectSize, abyHead, 4);
432 12 : CPL_MSBPTR32(&nSectSize);
433 12 : if (nSectSize < 5)
434 : {
435 0 : CPLDebug("GRIB", "Invalid section size for iSubMessage = %d",
436 : iSubMessage);
437 0 : return;
438 : }
439 12 : if (VSIFSeekL(poGDS->fp, nSectSize - 5, SEEK_CUR) != 0)
440 : {
441 0 : CPLDebug("GRIB",
442 : "Cannot read past section for iSubMessage = %d",
443 : iSubMessage);
444 0 : return;
445 : }
446 12 : if (abyHead[4] < 2 || abyHead[4] > 7)
447 : {
448 0 : CPLDebug("GRIB", "Invalid section number for iSubMessage = %d",
449 : iSubMessage);
450 0 : return;
451 : }
452 12 : if (abyHead[4] == 7)
453 : {
454 2 : ++iSubMessage;
455 : }
456 12 : if (VSIFReadL(abyHead, 5, 1, poGDS->fp) != 1)
457 : {
458 0 : CPLDebug("GRIB",
459 : "Cannot read 5 leading bytes for iSubMessage = %d",
460 : iSubMessage);
461 0 : return;
462 : }
463 : }
464 : }
465 :
466 : // Skip to section 4
467 930 : while (abyHead[4] != 4)
468 : {
469 611 : memcpy(&nSectSize, abyHead, 4);
470 611 : CPL_MSBPTR32(&nSectSize);
471 :
472 611 : const int nCurSection = abyHead[4];
473 611 : if (nSectSize < 5)
474 : {
475 0 : CPLDebug("GRIB", "Invalid section size for section %d",
476 : nCurSection);
477 0 : return;
478 : }
479 1222 : if (VSIFSeekL(poGDS->fp, nSectSize - 5, SEEK_CUR) != 0 ||
480 611 : VSIFReadL(abyHead, 5, 1, poGDS->fp) != 1)
481 : {
482 0 : CPLDebug("GRIB", "Cannot read section %d", nCurSection);
483 0 : return;
484 : }
485 : }
486 :
487 : // Collect section 4 octet information. We read the file
488 : // ourselves since the GRIB API does not appear to preserve all
489 : // this for us.
490 : // if( abyHead[4] == 4 )
491 : {
492 319 : memcpy(&nSectSize, abyHead, 4);
493 319 : CPL_MSBPTR32(&nSectSize);
494 319 : if (nSectSize >= 9 && nSectSize <= 100000 /* arbitrary upper limit */)
495 : {
496 319 : GByte *pabyBody = static_cast<GByte *>(CPLMalloc(nSectSize));
497 319 : memcpy(pabyBody, abyHead, 5);
498 319 : if (VSIFReadL(pabyBody + 5, 1, nSectSize - 5, poGDS->fp) !=
499 319 : nSectSize - 5)
500 : {
501 0 : CPLDebug("GRIB", "Cannot read section 4");
502 0 : CPLFree(pabyBody);
503 0 : return;
504 : }
505 :
506 319 : GUInt16 nCoordCount = 0;
507 319 : memcpy(&nCoordCount, pabyBody + 6 - 1, 2);
508 319 : CPL_MSBPTR16(&nCoordCount);
509 :
510 319 : GUInt16 nPDTN = 0;
511 319 : memcpy(&nPDTN, pabyBody + 8 - 1, 2);
512 319 : CPL_MSBPTR16(&nPDTN);
513 :
514 319 : GDALRasterBand::SetMetadataItem("GRIB_PDS_PDTN",
515 638 : CPLString().Printf("%d", nPDTN));
516 319 : m_nPDTN = nPDTN;
517 :
518 638 : CPLString osOctet;
519 319 : const int nTemplateFoundByteCount =
520 319 : nSectSize - 9U >= nCoordCount * 4U
521 319 : ? static_cast<int>(nSectSize - 9 - nCoordCount * 4)
522 : : 0;
523 9040 : for (int i = 0; i < nTemplateFoundByteCount; i++)
524 : {
525 8721 : char szByte[10] = {'\0'};
526 :
527 8721 : if (i == 0)
528 319 : snprintf(szByte, sizeof(szByte), "%d", pabyBody[i + 9]);
529 : else
530 8402 : snprintf(szByte, sizeof(szByte), " %d", pabyBody[i + 9]);
531 8721 : osOctet += szByte;
532 : }
533 :
534 319 : GDALRasterBand::SetMetadataItem("GRIB_PDS_TEMPLATE_NUMBERS",
535 : osOctet);
536 :
537 319 : g2int iofst = 0;
538 319 : g2int pdsnum = 0;
539 319 : g2int *pdstempl = nullptr;
540 319 : g2int mappdslen = 0;
541 319 : g2float *coordlist = nullptr;
542 319 : g2int numcoord = 0;
543 319 : if (getpdsindex(nPDTN) < 0)
544 : {
545 3 : CPLError(CE_Warning, CPLE_NotSupported,
546 : "Template 4.%d is not recognized currently", nPDTN);
547 : }
548 316 : else if (g2_unpack4(pabyBody, nSectSize, &iofst, &pdsnum, &pdstempl,
549 316 : &mappdslen, &coordlist, &numcoord) == 0)
550 : {
551 316 : gtemplate *mappds = extpdstemplate(pdsnum, pdstempl);
552 316 : if (mappds)
553 : {
554 316 : int nTemplateByteCount = 0;
555 5448 : for (int i = 0; i < mappds->maplen; i++)
556 5132 : nTemplateByteCount += abs(mappds->map[i]);
557 381 : for (int i = 0; i < mappds->extlen; i++)
558 65 : nTemplateByteCount += abs(mappds->ext[i]);
559 316 : if (nTemplateByteCount == nTemplateFoundByteCount)
560 : {
561 628 : CPLString osValues;
562 5479 : for (g2int i = 0; i < mappds->maplen + mappds->extlen;
563 : i++)
564 : {
565 5165 : if (i > 0)
566 4851 : osValues += " ";
567 5165 : const int nEltSize =
568 5165 : (i < mappds->maplen)
569 5165 : ? mappds->map[i]
570 65 : : mappds->ext[i - mappds->maplen];
571 5165 : if (nEltSize == 4)
572 : {
573 88 : m_anPDSTemplateAssembledValues.push_back(
574 88 : static_cast<GUInt32>(pdstempl[i]));
575 : osValues += CPLSPrintf(
576 88 : "%u", static_cast<GUInt32>(pdstempl[i]));
577 : }
578 : else
579 : {
580 5077 : m_anPDSTemplateAssembledValues.push_back(
581 5077 : pdstempl[i]);
582 5077 : osValues += CPLSPrintf("%d", pdstempl[i]);
583 : }
584 : }
585 314 : GDALRasterBand::SetMetadataItem(
586 : "GRIB_PDS_TEMPLATE_ASSEMBLED_VALUES", osValues);
587 : }
588 : else
589 : {
590 2 : CPLDebug(
591 : "GRIB",
592 : "Cannot expose GRIB_PDS_TEMPLATE_ASSEMBLED_VALUES "
593 : "as we would expect %d bytes from the "
594 : "tables, but %d are available",
595 : nTemplateByteCount, nTemplateFoundByteCount);
596 : }
597 :
598 316 : free(mappds->ext);
599 316 : free(mappds);
600 : }
601 : }
602 319 : free(pdstempl);
603 319 : free(coordlist);
604 :
605 319 : CPLFree(pabyBody);
606 :
607 638 : FindNoDataGrib2(false);
608 : }
609 : else
610 : {
611 0 : CPLDebug("GRIB", "Invalid section size for section %d", 4);
612 : }
613 : }
614 : }
615 :
616 : /************************************************************************/
617 : /* FindNoDataGrib2() */
618 : /************************************************************************/
619 :
620 319 : void GRIBRasterBand::FindNoDataGrib2(bool bSeekToStart)
621 : {
622 : // There is no easy way in the degrib API to retrieve the nodata value
623 : // without decompressing the data point section (which is slow), so
624 : // retrieve nodata value by parsing section 5 (Data Representation Section)
625 : // We also check section 6 to see if there is a bitmap
626 319 : GRIBDataset *poGDS = static_cast<GRIBDataset *>(poDS);
627 319 : CPLAssert(m_nGribVersion == 2);
628 :
629 319 : if (m_bHasLookedForNoData)
630 0 : return;
631 319 : m_bHasLookedForNoData = true;
632 :
633 319 : if (bSeekToStart)
634 : {
635 : // Skip over section 0
636 0 : VSIFSeekL(poGDS->fp, start + 16, SEEK_SET);
637 : }
638 :
639 319 : GByte abyHead[5] = {0};
640 319 : VSIFReadL(abyHead, 5, 1, poGDS->fp);
641 :
642 : // Skip to section 5
643 319 : GUInt32 nSectSize = 0;
644 319 : while (abyHead[4] != 5)
645 : {
646 0 : memcpy(&nSectSize, abyHead, 4);
647 0 : CPL_MSBPTR32(&nSectSize);
648 :
649 0 : if (nSectSize < 5 ||
650 0 : VSIFSeekL(poGDS->fp, nSectSize - 5, SEEK_CUR) != 0 ||
651 0 : VSIFReadL(abyHead, 5, 1, poGDS->fp) != 1)
652 0 : break;
653 : }
654 :
655 : // See http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_sect5.shtml
656 319 : if (abyHead[4] == 5)
657 : {
658 319 : memcpy(&nSectSize, abyHead, 4);
659 319 : CPL_MSBPTR32(&nSectSize);
660 319 : if (nSectSize >= 11 && nSectSize <= 100000 /* arbitrary upper limit */)
661 : {
662 319 : GByte *pabyBody = static_cast<GByte *>(CPLMalloc(nSectSize));
663 319 : memcpy(pabyBody, abyHead, 5);
664 319 : VSIFReadL(pabyBody + 5, 1, nSectSize - 5, poGDS->fp);
665 :
666 319 : GUInt16 nDRTN = 0;
667 319 : memcpy(&nDRTN, pabyBody + 10 - 1, 2);
668 319 : CPL_MSBPTR16(&nDRTN);
669 :
670 319 : GDALRasterBand::SetMetadataItem("DRS_DRTN", CPLSPrintf("%d", nDRTN),
671 : "GRIB");
672 319 : if ((nDRTN == GS5_SIMPLE || nDRTN == GS5_CMPLX ||
673 176 : nDRTN == GS5_CMPLXSEC || nDRTN == GS5_JPEG2000 ||
674 120 : nDRTN == GS5_PNG) &&
675 229 : nSectSize >= 20)
676 : {
677 : float fRef;
678 229 : memcpy(&fRef, pabyBody + 12 - 1, 4);
679 229 : CPL_MSBPTR32(&fRef);
680 229 : GDALRasterBand::SetMetadataItem(
681 : "DRS_REF_VALUE", CPLSPrintf("%.10f", fRef), "GRIB");
682 :
683 : GUInt16 nBinaryScaleFactorUnsigned;
684 229 : memcpy(&nBinaryScaleFactorUnsigned, pabyBody + 16 - 1, 2);
685 229 : CPL_MSBPTR16(&nBinaryScaleFactorUnsigned);
686 229 : const int nBSF =
687 : (nBinaryScaleFactorUnsigned & 0x8000)
688 25 : ? -static_cast<int>(nBinaryScaleFactorUnsigned & 0x7FFF)
689 229 : : static_cast<int>(nBinaryScaleFactorUnsigned);
690 229 : GDALRasterBand::SetMetadataItem("DRS_BINARY_SCALE_FACTOR",
691 : CPLSPrintf("%d", nBSF), "GRIB");
692 :
693 : GUInt16 nDecimalScaleFactorUnsigned;
694 229 : memcpy(&nDecimalScaleFactorUnsigned, pabyBody + 18 - 1, 2);
695 229 : CPL_MSBPTR16(&nDecimalScaleFactorUnsigned);
696 229 : const int nDSF =
697 : (nDecimalScaleFactorUnsigned & 0x8000)
698 19 : ? -static_cast<int>(nDecimalScaleFactorUnsigned &
699 : 0x7FFF)
700 229 : : static_cast<int>(nDecimalScaleFactorUnsigned);
701 229 : GDALRasterBand::SetMetadataItem("DRS_DECIMAL_SCALE_FACTOR",
702 : CPLSPrintf("%d", nDSF), "GRIB");
703 :
704 229 : const int nBits = pabyBody[20 - 1];
705 229 : GDALRasterBand::SetMetadataItem(
706 : "DRS_NBITS", CPLSPrintf("%d", nBits), "GRIB");
707 : }
708 :
709 : // 2 = Grid Point Data - Complex Packing
710 : // 3 = Grid Point Data - Complex Packing and Spatial Differencing
711 319 : if ((nDRTN == GS5_CMPLX || nDRTN == GS5_CMPLXSEC) &&
712 54 : nSectSize >= 31)
713 : {
714 54 : const int nMiss = pabyBody[23 - 1];
715 54 : if (nMiss == 1 || nMiss == 2)
716 : {
717 35 : const int original_field_type = pabyBody[21 - 1];
718 35 : if (original_field_type == 0) // Floating Point
719 : {
720 : float fTemp;
721 27 : memcpy(&fTemp, &pabyBody[24 - 1], 4);
722 27 : CPL_MSBPTR32(&fTemp);
723 27 : m_dfNoData = fTemp;
724 27 : m_bHasNoData = true;
725 27 : if (nMiss == 2)
726 : {
727 0 : memcpy(&fTemp, &pabyBody[28 - 1], 4);
728 0 : CPL_MSBPTR32(&fTemp);
729 0 : CPLDebug("GRIB",
730 : "Secondary missing value also set for "
731 : "band %d : %f",
732 : nBand, fTemp);
733 : }
734 : }
735 8 : else if (original_field_type == 1) // Integer
736 : {
737 : int iTemp;
738 8 : memcpy(&iTemp, &pabyBody[24 - 1], 4);
739 8 : CPL_MSBPTR32(&iTemp);
740 8 : m_dfNoData = iTemp;
741 8 : m_bHasNoData = true;
742 8 : if (nMiss == 2)
743 : {
744 0 : memcpy(&iTemp, &pabyBody[28 - 1], 4);
745 0 : CPL_MSBPTR32(&iTemp);
746 0 : CPLDebug("GRIB",
747 : "Secondary missing value also set for "
748 : "band %d : %d",
749 : nBand, iTemp);
750 : }
751 : }
752 : else
753 : {
754 : // FIXME What to do? Blindly convert to float?
755 0 : CPLDebug("GRIB",
756 : "Complex Packing - Type of Original Field "
757 : "Values for band %d: %u",
758 : nBand, original_field_type);
759 : }
760 : }
761 : }
762 :
763 319 : if (nDRTN == GS5_CMPLXSEC && nSectSize >= 48)
764 : {
765 21 : const int nOrder = pabyBody[48 - 1];
766 21 : GDALRasterBand::SetMetadataItem(
767 : "DRS_SPATIAL_DIFFERENCING_ORDER", CPLSPrintf("%d", nOrder),
768 : "GRIB");
769 : }
770 :
771 319 : CPLFree(pabyBody);
772 : }
773 0 : else if (nSectSize > 5)
774 : {
775 0 : VSIFSeekL(poGDS->fp, nSectSize - 5, SEEK_CUR);
776 : }
777 : }
778 :
779 319 : if (!m_bHasNoData)
780 : {
781 : // Check bitmap section
782 284 : GByte abySection6[6] = {0};
783 284 : VSIFReadL(abySection6, 6, 1, poGDS->fp);
784 : // Is there a bitmap ?
785 284 : if (abySection6[4] == 6 && abySection6[5] == 0)
786 : {
787 6 : m_dfNoData = 9999.0; // Same value as in metaparse.cpp:ParseGrid()
788 6 : m_bHasNoData = true;
789 : }
790 : }
791 : }
792 :
793 : /************************************************************************/
794 : /* GetDescription() */
795 : /************************************************************************/
796 :
797 23 : const char *GRIBRasterBand::GetDescription() const
798 : {
799 23 : if (longFstLevel == nullptr)
800 0 : return GDALPamRasterBand::GetDescription();
801 :
802 23 : return longFstLevel;
803 : }
804 :
805 : /************************************************************************/
806 : /* LoadData() */
807 : /************************************************************************/
808 :
809 9044 : CPLErr GRIBRasterBand::LoadData()
810 :
811 : {
812 9044 : if (!m_Grib_Data)
813 : {
814 178 : GRIBDataset *poGDS = static_cast<GRIBDataset *>(poDS);
815 :
816 178 : if (poGDS->bCacheOnlyOneBand)
817 : {
818 : // In "one-band-at-a-time" strategy, if the last recently used
819 : // band is not that one, uncache it. We could use a smarter strategy
820 : // based on a LRU, but that's a bit overkill for now.
821 0 : poGDS->poLastUsedBand->UncacheData();
822 0 : poGDS->nCachedBytes = 0;
823 : }
824 : else
825 : {
826 : // Once we have cached more than nCachedBytesThreshold bytes, we
827 : // will switch to "one-band-at-a-time" strategy, instead of caching
828 : // all bands that have been accessed.
829 178 : if (poGDS->nCachedBytes > poGDS->nCachedBytesThreshold)
830 : {
831 : GUIntBig nMinCacheSize =
832 0 : 1 + static_cast<GUIntBig>(poGDS->nRasterXSize) *
833 0 : poGDS->nRasterYSize * poGDS->nBands *
834 0 : GDALGetDataTypeSizeBytes(eDataType) / 1024 / 1024;
835 0 : CPLDebug("GRIB",
836 : "Maximum band cache size reached for this dataset. "
837 : "Caching only one band at a time from now, which can "
838 : "negatively affect performance. Consider "
839 : "increasing GRIB_CACHEMAX to a higher value (in MB), "
840 : "at least " CPL_FRMT_GUIB " in that instance",
841 : nMinCacheSize);
842 0 : for (int i = 0; i < poGDS->nBands; i++)
843 : {
844 : reinterpret_cast<GRIBRasterBand *>(
845 0 : poGDS->GetRasterBand(i + 1))
846 0 : ->UncacheData();
847 : }
848 0 : poGDS->nCachedBytes = 0;
849 0 : poGDS->bCacheOnlyOneBand = TRUE;
850 : }
851 : }
852 :
853 : // we don't seem to have any way to detect errors in this!
854 178 : if (m_Grib_MetaData != nullptr)
855 : {
856 133 : MetaFree(m_Grib_MetaData);
857 133 : delete m_Grib_MetaData;
858 133 : m_Grib_MetaData = nullptr;
859 : }
860 178 : ReadGribData(poGDS->fp, start, subgNum, &m_Grib_Data, &m_Grib_MetaData);
861 178 : if (!m_Grib_Data)
862 : {
863 0 : CPLError(CE_Failure, CPLE_AppDefined, "Out of memory.");
864 0 : if (m_Grib_MetaData != nullptr)
865 : {
866 0 : MetaFree(m_Grib_MetaData);
867 0 : delete m_Grib_MetaData;
868 0 : m_Grib_MetaData = nullptr;
869 : }
870 0 : return CE_Failure;
871 : }
872 :
873 : // Check the band matches the dataset as a whole, size wise. (#3246)
874 178 : nGribDataXSize = m_Grib_MetaData->gds.Nx;
875 178 : nGribDataYSize = m_Grib_MetaData->gds.Ny;
876 178 : if (nGribDataXSize <= 0 || nGribDataYSize <= 0)
877 : {
878 0 : CPLError(CE_Failure, CPLE_AppDefined,
879 : "Band %d of GRIB dataset is %dx%d.", nBand, nGribDataXSize,
880 : nGribDataYSize);
881 0 : MetaFree(m_Grib_MetaData);
882 0 : delete m_Grib_MetaData;
883 0 : m_Grib_MetaData = nullptr;
884 0 : return CE_Failure;
885 : }
886 :
887 178 : poGDS->nCachedBytes += static_cast<GIntBig>(nGribDataXSize) *
888 178 : nGribDataYSize * sizeof(double);
889 178 : poGDS->poLastUsedBand = this;
890 :
891 178 : if (nGribDataXSize != nRasterXSize || nGribDataYSize != nRasterYSize)
892 : {
893 24 : CPLError(CE_Warning, CPLE_AppDefined,
894 : "Band %d of GRIB dataset is %dx%d, while the first band "
895 : "and dataset is %dx%d. Georeferencing of band %d may "
896 : "be incorrect, and data access may be incomplete.",
897 : nBand, nGribDataXSize, nGribDataYSize, nRasterXSize,
898 : nRasterYSize, nBand);
899 : }
900 : }
901 :
902 9044 : return CE_None;
903 : }
904 :
905 : /************************************************************************/
906 : /* GetMetaData() */
907 : /************************************************************************/
908 86 : char **GRIBRasterBand::GetMetadata(const char *pszDomain)
909 : {
910 86 : FindMetaData();
911 138 : if (m_nGribVersion == 2 &&
912 52 : CPLTestBool(CPLGetConfigOption("GRIB_PDS_ALL_BANDS", "ON")))
913 : {
914 52 : FindPDSTemplateGRIB2();
915 : }
916 86 : return GDALPamRasterBand::GetMetadata(pszDomain);
917 : }
918 :
919 : /************************************************************************/
920 : /* GetMetaDataItem() */
921 : /************************************************************************/
922 693 : const char *GRIBRasterBand::GetMetadataItem(const char *pszName,
923 : const char *pszDomain)
924 : {
925 693 : FindMetaData();
926 1231 : if (m_nGribVersion == 2 &&
927 538 : CPLTestBool(CPLGetConfigOption("GRIB_PDS_ALL_BANDS", "ON")))
928 : {
929 537 : FindPDSTemplateGRIB2();
930 : }
931 693 : return GDALPamRasterBand::GetMetadataItem(pszName, pszDomain);
932 : }
933 :
934 : /************************************************************************/
935 : /* IReadBlock() */
936 : /************************************************************************/
937 :
938 9044 : CPLErr GRIBRasterBand::IReadBlock(int /* nBlockXOff */, int nBlockYOff,
939 : void *pImage)
940 :
941 : {
942 9044 : CPLErr eErr = LoadData();
943 9044 : if (eErr != CE_None)
944 0 : return eErr;
945 :
946 9044 : GRIBDataset *poGDS = static_cast<GRIBDataset *>(poDS);
947 :
948 : // The image as read is always upside down to our normal
949 : // orientation so we need to effectively flip it at this
950 : // point. We also need to deal with bands that are a different
951 : // size than the dataset as a whole.
952 :
953 9044 : if (nGribDataXSize == nRasterXSize && nGribDataYSize == nRasterYSize &&
954 8552 : poGDS->nSplitAndSwapColumn == 0)
955 : {
956 : // Simple 1:1 case.
957 6893 : memcpy(pImage,
958 6893 : m_Grib_Data + static_cast<size_t>(nRasterXSize) *
959 6893 : (nRasterYSize - nBlockYOff - 1),
960 6893 : nRasterXSize * sizeof(double));
961 :
962 6893 : return CE_None;
963 : }
964 :
965 2151 : memset(pImage, 0, sizeof(double) * nRasterXSize);
966 :
967 2151 : if (nBlockYOff >= nGribDataYSize) // Off image?
968 57 : return CE_None;
969 :
970 2094 : int nSplitAndSwapColumn = poGDS->nSplitAndSwapColumn;
971 2094 : if (nRasterXSize != nGribDataXSize)
972 435 : nSplitAndSwapColumn = 0;
973 :
974 2094 : const int nCopyWords = std::min(nRasterXSize, nGribDataXSize);
975 :
976 2094 : memcpy(pImage,
977 2094 : m_Grib_Data +
978 2094 : static_cast<size_t>(nGribDataXSize) *
979 2094 : (nGribDataYSize - nBlockYOff - 1) +
980 : nSplitAndSwapColumn,
981 2094 : (nCopyWords - nSplitAndSwapColumn) * sizeof(double));
982 :
983 2094 : if (nSplitAndSwapColumn > 0)
984 1659 : memcpy(reinterpret_cast<void *>(reinterpret_cast<double *>(pImage) +
985 1659 : nCopyWords - nSplitAndSwapColumn),
986 1659 : m_Grib_Data + static_cast<size_t>(nGribDataXSize) *
987 1659 : (nGribDataYSize - nBlockYOff - 1),
988 1659 : nSplitAndSwapColumn * sizeof(double));
989 :
990 2094 : return CE_None;
991 : }
992 :
993 : /************************************************************************/
994 : /* GetNoDataValue() */
995 : /************************************************************************/
996 :
997 125 : double GRIBRasterBand::GetNoDataValue(int *pbSuccess)
998 : {
999 125 : if (m_bHasLookedForNoData)
1000 : {
1001 105 : if (pbSuccess)
1002 105 : *pbSuccess = m_bHasNoData;
1003 105 : return m_dfNoData;
1004 : }
1005 :
1006 20 : m_bHasLookedForNoData = true;
1007 20 : if (m_Grib_MetaData == nullptr)
1008 : {
1009 18 : GRIBDataset *poGDS = static_cast<GRIBDataset *>(poDS);
1010 18 : ReadGribData(poGDS->fp, start, subgNum, nullptr, &m_Grib_MetaData);
1011 18 : if (m_Grib_MetaData == nullptr)
1012 : {
1013 0 : m_bHasNoData = false;
1014 0 : m_dfNoData = 0;
1015 0 : if (pbSuccess)
1016 0 : *pbSuccess = m_bHasNoData;
1017 0 : return m_dfNoData;
1018 : }
1019 : }
1020 :
1021 20 : if (m_Grib_MetaData->gridAttrib.f_miss == 0)
1022 : {
1023 4 : m_bHasNoData = false;
1024 4 : m_dfNoData = 0;
1025 4 : if (pbSuccess)
1026 4 : *pbSuccess = m_bHasNoData;
1027 4 : return m_dfNoData;
1028 : }
1029 :
1030 16 : if (m_Grib_MetaData->gridAttrib.f_miss == 2)
1031 : {
1032 : // What TODO?
1033 0 : CPLDebug("GRIB", "Secondary missing value also set for band %d : %f",
1034 0 : nBand, m_Grib_MetaData->gridAttrib.missSec);
1035 : }
1036 :
1037 16 : m_bHasNoData = true;
1038 16 : m_dfNoData = m_Grib_MetaData->gridAttrib.missPri;
1039 16 : if (pbSuccess)
1040 16 : *pbSuccess = m_bHasNoData;
1041 16 : return m_dfNoData;
1042 : }
1043 :
1044 : /************************************************************************/
1045 : /* ReadGribData() */
1046 : /************************************************************************/
1047 :
1048 523 : void GRIBRasterBand::ReadGribData(VSILFILE *fp, vsi_l_offset start, int subgNum,
1049 : double **data, grib_MetaData **metaData)
1050 : {
1051 : // Initialization, for calling the ReadGrib2Record function.
1052 523 : sInt4 f_endMsg = 1; // 1 if we read the last grid in a GRIB message, or we
1053 : // haven't read any messages.
1054 : // int subgNum = 0; // The subgrid in the message that we are interested in.
1055 523 : sChar f_unit = 2; // None = 0, English = 1, Metric = 2
1056 523 : double majEarth = 0.0; // -radEarth if < 6000 ignore, otherwise use this
1057 : // to override the radEarth in the GRIB1 or GRIB2
1058 : // message. Needed because NCEP uses 6371.2 but
1059 : // GRIB1 could only state 6367.47.
1060 523 : double minEarth = 0.0; // -minEarth if < 6000 ignore, otherwise use this
1061 : // to override the minEarth in the GRIB1 or GRIB2
1062 : // message.
1063 523 : sChar f_SimpleVer = 4; // Which version of the simple NDFD Weather table
1064 : // to use. (1 is 6/2003) (2 is 1/2004) (3 is
1065 : // 2/2004) (4 is 11/2004) (default 4)
1066 : LatLon lwlf; // Lower left corner (cookie slicing) -lwlf
1067 : LatLon uprt; // Upper right corner (cookie slicing) -uprt
1068 : IS_dataType is; // Un-parsed meta data for this GRIB2 message. As well
1069 : // as some memory used by the unpacker.
1070 :
1071 523 : lwlf.lat = -100; // lat == -100 instructs the GRIB decoder that we don't
1072 : // want a subgrid
1073 :
1074 523 : IS_Init(&is);
1075 :
1076 : const char *pszGribNormalizeUnits =
1077 523 : CPLGetConfigOption("GRIB_NORMALIZE_UNITS", "YES");
1078 523 : if (!CPLTestBool(pszGribNormalizeUnits))
1079 2 : f_unit = 0; // Do not normalize units to metric.
1080 :
1081 523 : start = FindTrueStart(fp, start);
1082 : // Read GRIB message from file position "start".
1083 523 : VSIFSeekL(fp, start, SEEK_SET);
1084 523 : uInt4 grib_DataLen = 0; // Size of Grib_Data.
1085 523 : *metaData = new grib_MetaData();
1086 523 : MetaInit(*metaData);
1087 523 : const int simpWWA = 0; // seem to be unused in degrib
1088 523 : ReadGrib2Record(fp, f_unit, data, &grib_DataLen, *metaData, &is, subgNum,
1089 : majEarth, minEarth, f_SimpleVer, simpWWA, &f_endMsg, &lwlf,
1090 : &uprt);
1091 :
1092 : // No intention to show errors, just swallow it and free the memory.
1093 523 : char *errMsg = errSprintf(nullptr);
1094 523 : if (errMsg != nullptr)
1095 339 : CPLDebug("GRIB", "%s", errMsg);
1096 523 : free(errMsg);
1097 523 : IS_Free(&is);
1098 523 : }
1099 :
1100 : /************************************************************************/
1101 : /* UncacheData() */
1102 : /************************************************************************/
1103 :
1104 431 : void GRIBRasterBand::UncacheData()
1105 : {
1106 431 : if (m_Grib_Data)
1107 178 : free(m_Grib_Data);
1108 431 : m_Grib_Data = nullptr;
1109 431 : if (m_Grib_MetaData)
1110 : {
1111 362 : MetaFree(m_Grib_MetaData);
1112 362 : delete m_Grib_MetaData;
1113 : }
1114 431 : m_Grib_MetaData = nullptr;
1115 431 : }
1116 :
1117 : /************************************************************************/
1118 : /* ~GRIBRasterBand() */
1119 : /************************************************************************/
1120 :
1121 839 : GRIBRasterBand::~GRIBRasterBand()
1122 : {
1123 431 : if (longFstLevel != nullptr)
1124 431 : CPLFree(longFstLevel);
1125 431 : UncacheData();
1126 839 : }
1127 :
1128 : /************************************************************************/
1129 : /* InventoryWrapperGrib */
1130 : /************************************************************************/
1131 : class InventoryWrapperGrib : public gdal::grib::InventoryWrapper
1132 : {
1133 : public:
1134 304 : explicit InventoryWrapperGrib(VSILFILE *fp) : gdal::grib::InventoryWrapper()
1135 : {
1136 304 : result_ = GRIB2Inventory(fp, &inv_, &inv_len_, 0 /* all messages */,
1137 : &num_messages_);
1138 304 : }
1139 :
1140 608 : ~InventoryWrapperGrib() override
1141 304 : {
1142 304 : if (inv_ == nullptr)
1143 0 : return;
1144 719 : for (uInt4 i = 0; i < inv_len_; i++)
1145 : {
1146 415 : GRIB2InventoryFree(inv_ + i);
1147 : }
1148 304 : free(inv_);
1149 608 : }
1150 : };
1151 :
1152 : /************************************************************************/
1153 : /* InventoryWrapperSidecar */
1154 : /************************************************************************/
1155 :
1156 : class InventoryWrapperSidecar : public gdal::grib::InventoryWrapper
1157 : {
1158 : public:
1159 3 : explicit InventoryWrapperSidecar(VSILFILE *fp)
1160 3 : : gdal::grib::InventoryWrapper()
1161 : {
1162 3 : result_ = -1;
1163 3 : VSIFSeekL(fp, 0, SEEK_END);
1164 3 : size_t length = static_cast<size_t>(VSIFTellL(fp));
1165 3 : if (length > 4 * 1024 * 1024)
1166 0 : return;
1167 3 : std::string psSidecar;
1168 3 : psSidecar.resize(length);
1169 3 : VSIFSeekL(fp, 0, SEEK_SET);
1170 3 : if (VSIFReadL(&psSidecar[0], length, 1, fp) != 1)
1171 0 : return;
1172 :
1173 : CPLStringList aosMsgs(
1174 : CSLTokenizeString2(psSidecar.c_str(), "\n",
1175 3 : CSLT_PRESERVEQUOTES | CSLT_STRIPLEADSPACES));
1176 3 : inv_len_ = aosMsgs.size();
1177 3 : inv_ = static_cast<inventoryType *>(
1178 3 : CPLMalloc(inv_len_ * sizeof(inventoryType)));
1179 :
1180 29 : for (size_t i = 0; i < inv_len_; ++i)
1181 : {
1182 : // We are parsing
1183 : // "msgNum[.subgNum]:start:dontcare:name1:name2:name3" For NOMADS:
1184 : // "msgNum[.subgNum]:start:reftime:var:level:time"
1185 : CPLStringList aosTokens(CSLTokenizeString2(
1186 26 : aosMsgs[i], ":", CSLT_PRESERVEQUOTES | CSLT_ALLOWEMPTYTOKENS));
1187 26 : CPLStringList aosNum;
1188 :
1189 26 : if (aosTokens.size() < 6)
1190 0 : goto err_sidecar;
1191 :
1192 26 : aosNum = CPLStringList(CSLTokenizeString2(aosTokens[0], ".", 0));
1193 26 : if (aosNum.size() < 1)
1194 0 : goto err_sidecar;
1195 :
1196 : // FindMetaData will retrieve the correct version number
1197 : char *endptr;
1198 26 : strtol(aosNum[0], &endptr, 10);
1199 26 : if (*endptr != 0)
1200 0 : goto err_sidecar;
1201 :
1202 26 : if (aosNum.size() < 2)
1203 24 : inv_[i].subgNum = 0;
1204 : else
1205 : {
1206 2 : auto subgNum = strtol(aosNum[1], &endptr, 10);
1207 2 : if (*endptr != 0)
1208 0 : goto err_sidecar;
1209 2 : if (subgNum <= 0 || subgNum > 65536)
1210 0 : goto err_sidecar;
1211 : // .idx file use a 1-based indexing, whereas DEGRIB uses a
1212 : // 0-based one
1213 2 : subgNum--;
1214 2 : inv_[i].subgNum = static_cast<unsigned short>(subgNum);
1215 : }
1216 :
1217 26 : inv_[i].start = strtoll(aosTokens[1], &endptr, 10);
1218 26 : if (*endptr != 0)
1219 0 : goto err_sidecar;
1220 :
1221 26 : inv_[i].unitName = nullptr;
1222 26 : inv_[i].comment = nullptr;
1223 26 : inv_[i].element = nullptr;
1224 26 : inv_[i].shortFstLevel = nullptr;
1225 : // This is going into the description field ->
1226 : // the only one available before loading the metadata
1227 26 : inv_[i].longFstLevel = VSIStrdup(CPLSPrintf(
1228 : "%s:%s:%s", aosTokens[3], aosTokens[4], aosTokens[5]));
1229 :
1230 26 : continue;
1231 :
1232 0 : err_sidecar:
1233 0 : CPLDebug("GRIB",
1234 : "Failed parsing sidecar entry '%s', "
1235 : "falling back to constructing an inventory",
1236 : aosMsgs[i]);
1237 0 : inv_len_ = static_cast<unsigned>(i);
1238 0 : return;
1239 : }
1240 :
1241 3 : result_ = inv_len_;
1242 : }
1243 :
1244 6 : ~InventoryWrapperSidecar() override
1245 3 : {
1246 3 : if (inv_ == nullptr)
1247 0 : return;
1248 :
1249 29 : for (unsigned i = 0; i < inv_len_; i++)
1250 26 : VSIFree(inv_[i].longFstLevel);
1251 :
1252 3 : VSIFree(inv_);
1253 6 : }
1254 : };
1255 :
1256 : /************************************************************************/
1257 : /* ==================================================================== */
1258 : /* GRIBDataset */
1259 : /* ==================================================================== */
1260 : /************************************************************************/
1261 :
1262 307 : GRIBDataset::GRIBDataset()
1263 : : fp(nullptr), nCachedBytes(0),
1264 : // Switch caching strategy once 100 MB threshold is reached.
1265 : // Why 100 MB? --> Why not.
1266 614 : nCachedBytesThreshold(static_cast<GIntBig>(atoi(
1267 : CPLGetConfigOption("GRIB_CACHEMAX", "100"))) *
1268 307 : 1024 * 1024),
1269 307 : bCacheOnlyOneBand(FALSE), nSplitAndSwapColumn(0), poLastUsedBand(nullptr)
1270 : {
1271 307 : adfGeoTransform[0] = 0.0;
1272 307 : adfGeoTransform[1] = 1.0;
1273 307 : adfGeoTransform[2] = 0.0;
1274 307 : adfGeoTransform[3] = 0.0;
1275 307 : adfGeoTransform[4] = 0.0;
1276 307 : adfGeoTransform[5] = 1.0;
1277 307 : }
1278 :
1279 : /************************************************************************/
1280 : /* ~GRIBDataset() */
1281 : /************************************************************************/
1282 :
1283 614 : GRIBDataset::~GRIBDataset()
1284 :
1285 : {
1286 307 : FlushCache(true);
1287 307 : if (fp != nullptr)
1288 303 : VSIFCloseL(fp);
1289 614 : }
1290 :
1291 : /************************************************************************/
1292 : /* GetGeoTransform() */
1293 : /************************************************************************/
1294 :
1295 133 : CPLErr GRIBDataset::GetGeoTransform(double *padfTransform)
1296 :
1297 : {
1298 133 : memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
1299 133 : return CE_None;
1300 : }
1301 :
1302 : /************************************************************************/
1303 : /* Inventory() */
1304 : /************************************************************************/
1305 :
1306 : std::unique_ptr<gdal::grib::InventoryWrapper>
1307 303 : GRIBDataset::Inventory(VSILFILE *fp, GDALOpenInfo *poOpenInfo)
1308 : {
1309 303 : std::unique_ptr<gdal::grib::InventoryWrapper> pInventories;
1310 :
1311 303 : VSIFSeekL(fp, 0, SEEK_SET);
1312 909 : CPLString sSideCarFilename = CPLString(poOpenInfo->pszFilename) + ".idx";
1313 303 : VSILFILE *fpSideCar = nullptr;
1314 303 : if (CPLTestBool(CSLFetchNameValueDef(poOpenInfo->papszOpenOptions,
1315 603 : "USE_IDX", "YES")) &&
1316 300 : ((fpSideCar = VSIFOpenL(sSideCarFilename, "rb")) != nullptr))
1317 : {
1318 3 : CPLDebug("GRIB", "Reading inventories from sidecar file %s",
1319 : sSideCarFilename.c_str());
1320 : // Contains an GRIB2 message inventory of the file.
1321 3 : pInventories = std::make_unique<InventoryWrapperSidecar>(fpSideCar);
1322 3 : if (pInventories->result() <= 0 || pInventories->length() == 0)
1323 0 : pInventories = nullptr;
1324 3 : VSIFCloseL(fpSideCar);
1325 : }
1326 : else
1327 300 : CPLDebug("GRIB", "Failed opening sidecar %s", sSideCarFilename.c_str());
1328 :
1329 303 : if (pInventories == nullptr)
1330 : {
1331 300 : CPLDebug("GRIB", "Reading inventories from GRIB file %s",
1332 : poOpenInfo->pszFilename);
1333 : // Contains an GRIB2 message inventory of the file.
1334 300 : pInventories = std::make_unique<InventoryWrapperGrib>(fp);
1335 : }
1336 :
1337 606 : return pInventories;
1338 : }
1339 :
1340 : /************************************************************************/
1341 : /* Open() */
1342 : /************************************************************************/
1343 :
1344 309 : GDALDataset *GRIBDataset::Open(GDALOpenInfo *poOpenInfo)
1345 :
1346 : {
1347 : #ifndef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
1348 : // During fuzzing, do not use Identify to reject crazy content.
1349 309 : if (!GRIBDriverIdentify(poOpenInfo))
1350 2 : return nullptr;
1351 : #endif
1352 307 : if (poOpenInfo->fpL == nullptr)
1353 0 : return nullptr;
1354 :
1355 : // A fast "probe" on the header that is partially read in memory.
1356 307 : char *buff = nullptr;
1357 307 : uInt4 buffLen = 0;
1358 307 : sInt4 sect0[SECT0LEN_WORD] = {0};
1359 307 : uInt4 gribLen = 0;
1360 307 : int version = 0;
1361 :
1362 : // grib is not thread safe, make sure not to cause problems
1363 : // for other thread safe formats
1364 614 : CPLMutexHolderD(&hGRIBMutex);
1365 :
1366 614 : CPLString tmpFilename;
1367 307 : tmpFilename.Printf("/vsimem/gribdataset-%p", poOpenInfo);
1368 :
1369 307 : VSILFILE *memfp = VSIFileFromMemBuffer(tmpFilename, poOpenInfo->pabyHeader,
1370 307 : poOpenInfo->nHeaderBytes, FALSE);
1371 614 : if (memfp == nullptr ||
1372 307 : ReadSECT0(memfp, &buff, &buffLen, -1, sect0, &gribLen, &version) < 0)
1373 : {
1374 0 : if (memfp != nullptr)
1375 : {
1376 0 : VSIFCloseL(memfp);
1377 0 : VSIUnlink(tmpFilename);
1378 : }
1379 0 : free(buff);
1380 0 : char *errMsg = errSprintf(nullptr);
1381 0 : if (errMsg != nullptr && strstr(errMsg, "Ran out of file") == nullptr)
1382 0 : CPLDebug("GRIB", "%s", errMsg);
1383 0 : free(errMsg);
1384 0 : return nullptr;
1385 : }
1386 307 : VSIFCloseL(memfp);
1387 307 : VSIUnlink(tmpFilename);
1388 307 : free(buff);
1389 :
1390 : // Confirm the requested access is supported.
1391 307 : if (poOpenInfo->eAccess == GA_Update)
1392 : {
1393 0 : CPLError(CE_Failure, CPLE_NotSupported,
1394 : "The GRIB driver does not support update access to existing "
1395 : "datasets.");
1396 0 : return nullptr;
1397 : }
1398 :
1399 307 : if (poOpenInfo->nOpenFlags & GDAL_OF_MULTIDIM_RASTER)
1400 : {
1401 4 : return OpenMultiDim(poOpenInfo);
1402 : }
1403 :
1404 : // Create a corresponding GDALDataset.
1405 303 : GRIBDataset *poDS = new GRIBDataset();
1406 :
1407 303 : poDS->fp = poOpenInfo->fpL;
1408 303 : poOpenInfo->fpL = nullptr;
1409 :
1410 : // Make an inventory of the GRIB file.
1411 : // The inventory does not contain all the information needed for
1412 : // creating the RasterBands (especially the x and y size), therefore
1413 : // the first GRIB band is also read for some additional metadata.
1414 : // The band-data that is read is stored into the first RasterBand,
1415 : // simply so that the same portion of the file is not read twice.
1416 :
1417 606 : auto pInventories = Inventory(poDS->fp, poOpenInfo);
1418 303 : if (pInventories->result() <= 0)
1419 : {
1420 9 : char *errMsg = errSprintf(nullptr);
1421 9 : if (errMsg != nullptr)
1422 8 : CPLDebug("GRIB", "%s", errMsg);
1423 9 : free(errMsg);
1424 :
1425 9 : CPLError(CE_Failure, CPLE_OpenFailed,
1426 : "%s is a grib file, "
1427 : "but no raster dataset was successfully identified.",
1428 : poOpenInfo->pszFilename);
1429 : // Release hGRIBMutex otherwise we'll deadlock with GDALDataset own
1430 : // hGRIBMutex.
1431 9 : CPLReleaseMutex(hGRIBMutex);
1432 9 : delete poDS;
1433 9 : CPLAcquireMutex(hGRIBMutex, 1000.0);
1434 9 : return nullptr;
1435 : }
1436 :
1437 : // Create band objects.
1438 702 : for (uInt4 i = 0; i < pInventories->length(); ++i)
1439 : {
1440 408 : inventoryType *psInv = pInventories->get(i);
1441 408 : GRIBRasterBand *gribBand = nullptr;
1442 408 : uInt4 bandNr = i + 1;
1443 :
1444 408 : if (bandNr == 1)
1445 : {
1446 : // Important: set DataSet extents before creating first RasterBand
1447 : // in it.
1448 294 : grib_MetaData *metaData = nullptr;
1449 294 : GRIBRasterBand::ReadGribData(poDS->fp, 0, psInv->subgNum, nullptr,
1450 : &metaData);
1451 294 : if (metaData == nullptr || metaData->gds.Nx < 1 ||
1452 294 : metaData->gds.Ny < 1)
1453 : {
1454 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1455 : "%s is a grib file, "
1456 : "but no raster dataset was successfully identified.",
1457 : poOpenInfo->pszFilename);
1458 : // Release hGRIBMutex otherwise we'll deadlock with GDALDataset
1459 : // own hGRIBMutex.
1460 0 : CPLReleaseMutex(hGRIBMutex);
1461 0 : delete poDS;
1462 0 : CPLAcquireMutex(hGRIBMutex, 1000.0);
1463 0 : if (metaData != nullptr)
1464 : {
1465 0 : MetaFree(metaData);
1466 0 : delete metaData;
1467 : }
1468 0 : return nullptr;
1469 : }
1470 294 : psInv->GribVersion = metaData->GribVersion;
1471 :
1472 : // Set the DataSet's x,y size, georeference and projection from
1473 : // the first GRIB band.
1474 294 : poDS->SetGribMetaData(metaData);
1475 294 : gribBand = new GRIBRasterBand(poDS, bandNr, psInv);
1476 :
1477 294 : if (psInv->GribVersion == 2)
1478 287 : gribBand->FindPDSTemplateGRIB2();
1479 :
1480 294 : gribBand->m_Grib_MetaData = metaData;
1481 : }
1482 : else
1483 : {
1484 114 : gribBand = new GRIBRasterBand(poDS, bandNr, psInv);
1485 : }
1486 408 : poDS->SetBand(bandNr, gribBand);
1487 : }
1488 :
1489 : // Initialize any PAM information.
1490 294 : poDS->SetDescription(poOpenInfo->pszFilename);
1491 :
1492 : // Release hGRIBMutex otherwise we'll deadlock with GDALDataset own
1493 : // hGRIBMutex.
1494 294 : CPLReleaseMutex(hGRIBMutex);
1495 294 : poDS->TryLoadXML();
1496 :
1497 : // Check for external overviews.
1498 294 : poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename,
1499 : poOpenInfo->GetSiblingFiles());
1500 294 : CPLAcquireMutex(hGRIBMutex, 1000.0);
1501 :
1502 294 : return poDS;
1503 : }
1504 :
1505 : /************************************************************************/
1506 : /* GRIBSharedResource */
1507 : /************************************************************************/
1508 :
1509 : struct GRIBSharedResource
1510 : {
1511 : VSILFILE *m_fp = nullptr;
1512 : vsi_l_offset m_nOffsetCurData = static_cast<vsi_l_offset>(-1);
1513 : std::vector<double> m_adfCurData{};
1514 : std::string m_osFilename;
1515 : std::shared_ptr<GDALPamMultiDim> m_poPAM{};
1516 :
1517 : GRIBSharedResource(const std::string &osFilename, VSILFILE *fp);
1518 : ~GRIBSharedResource();
1519 :
1520 : const std::vector<double> &LoadData(vsi_l_offset nOffset, int subgNum);
1521 :
1522 23 : const std::shared_ptr<GDALPamMultiDim> &GetPAM()
1523 : {
1524 23 : return m_poPAM;
1525 : }
1526 : };
1527 :
1528 4 : GRIBSharedResource::GRIBSharedResource(const std::string &osFilename,
1529 4 : VSILFILE *fp)
1530 : : m_fp(fp), m_osFilename(osFilename),
1531 4 : m_poPAM(std::make_shared<GDALPamMultiDim>(osFilename))
1532 : {
1533 4 : }
1534 :
1535 4 : GRIBSharedResource::~GRIBSharedResource()
1536 : {
1537 4 : if (m_fp)
1538 4 : VSIFCloseL(m_fp);
1539 4 : }
1540 :
1541 : /************************************************************************/
1542 : /* GRIBGroup */
1543 : /************************************************************************/
1544 :
1545 : class GRIBArray;
1546 :
1547 : class GRIBGroup final : public GDALGroup
1548 : {
1549 : friend class GRIBArray;
1550 : std::shared_ptr<GRIBSharedResource> m_poShared{};
1551 : std::vector<std::shared_ptr<GDALMDArray>> m_poArrays{};
1552 : std::vector<std::shared_ptr<GDALDimension>> m_dims{};
1553 : std::map<std::string, std::shared_ptr<GDALDimension>> m_oMapDims{};
1554 : int m_nHorizDimCounter = 0;
1555 : std::shared_ptr<GDALGroup> m_memRootGroup{};
1556 :
1557 : public:
1558 4 : explicit GRIBGroup(const std::shared_ptr<GRIBSharedResource> &poShared)
1559 4 : : GDALGroup(std::string(), "/"), m_poShared(poShared)
1560 : {
1561 : std::unique_ptr<GDALDataset> poTmpDS(
1562 4 : MEMDataset::CreateMultiDimensional("", nullptr, nullptr));
1563 4 : m_memRootGroup = poTmpDS->GetRootGroup();
1564 4 : }
1565 :
1566 36 : void AddArray(const std::shared_ptr<GDALMDArray> &array)
1567 : {
1568 36 : m_poArrays.emplace_back(array);
1569 36 : }
1570 :
1571 : std::vector<std::string>
1572 : GetMDArrayNames(CSLConstList papszOptions) const override;
1573 : std::shared_ptr<GDALMDArray>
1574 : OpenMDArray(const std::string &osName,
1575 : CSLConstList papszOptions) const override;
1576 :
1577 : std::vector<std::shared_ptr<GDALDimension>>
1578 3 : GetDimensions(CSLConstList) const override
1579 : {
1580 3 : return m_dims;
1581 : }
1582 : };
1583 :
1584 : /************************************************************************/
1585 : /* GRIBArray */
1586 : /************************************************************************/
1587 :
1588 : class GRIBArray final : public GDALPamMDArray
1589 : {
1590 : std::shared_ptr<GRIBSharedResource> m_poShared;
1591 : std::vector<std::shared_ptr<GDALDimension>> m_dims{};
1592 : GDALExtendedDataType m_dt = GDALExtendedDataType::Create(GDT_Float64);
1593 : std::shared_ptr<OGRSpatialReference> m_poSRS{};
1594 : std::vector<vsi_l_offset> m_anOffsets{};
1595 : std::vector<int> m_anSubgNums{};
1596 : std::vector<double> m_adfTimes{};
1597 : std::vector<std::shared_ptr<GDALAttribute>> m_attributes{};
1598 : std::string m_osUnit{};
1599 : std::vector<GByte> m_abyNoData{};
1600 :
1601 : GRIBArray(const std::string &osName,
1602 : const std::shared_ptr<GRIBSharedResource> &poShared);
1603 :
1604 : protected:
1605 : bool IRead(const GUInt64 *arrayStartIdx, const size_t *count,
1606 : const GInt64 *arrayStep, const GPtrDiff_t *bufferStride,
1607 : const GDALExtendedDataType &bufferDataType,
1608 : void *pDstBuffer) const override;
1609 :
1610 : public:
1611 : static std::shared_ptr<GRIBArray>
1612 23 : Create(const std::string &osName,
1613 : const std::shared_ptr<GRIBSharedResource> &poShared)
1614 : {
1615 23 : auto ar(std::shared_ptr<GRIBArray>(new GRIBArray(osName, poShared)));
1616 23 : ar->SetSelf(ar);
1617 23 : return ar;
1618 : }
1619 :
1620 : void Init(GRIBGroup *poGroup, GRIBDataset *poDS, GRIBRasterBand *poBand,
1621 : inventoryType *psInv);
1622 : void ExtendTimeDim(vsi_l_offset nOffset, int subgNum, double dfValidTime);
1623 : void Finalize(GRIBGroup *poGroup, inventoryType *psInv);
1624 :
1625 0 : bool IsWritable() const override
1626 : {
1627 0 : return false;
1628 : }
1629 :
1630 4 : const std::string &GetFilename() const override
1631 : {
1632 4 : return m_poShared->m_osFilename;
1633 : }
1634 :
1635 : const std::vector<std::shared_ptr<GDALDimension>> &
1636 30 : GetDimensions() const override
1637 : {
1638 30 : return m_dims;
1639 : }
1640 :
1641 11 : const GDALExtendedDataType &GetDataType() const override
1642 : {
1643 11 : return m_dt;
1644 : }
1645 :
1646 1 : std::shared_ptr<OGRSpatialReference> GetSpatialRef() const override
1647 : {
1648 1 : return m_poSRS;
1649 : }
1650 :
1651 : std::vector<std::shared_ptr<GDALAttribute>>
1652 14 : GetAttributes(CSLConstList) const override
1653 : {
1654 14 : return m_attributes;
1655 : }
1656 :
1657 1 : const std::string &GetUnit() const override
1658 : {
1659 1 : return m_osUnit;
1660 : }
1661 :
1662 1 : const void *GetRawNoDataValue() const override
1663 : {
1664 1 : return m_abyNoData.empty() ? nullptr : m_abyNoData.data();
1665 : }
1666 : };
1667 :
1668 : /************************************************************************/
1669 : /* GetMDArrayNames() */
1670 : /************************************************************************/
1671 :
1672 3 : std::vector<std::string> GRIBGroup::GetMDArrayNames(CSLConstList) const
1673 : {
1674 3 : std::vector<std::string> ret;
1675 31 : for (const auto &array : m_poArrays)
1676 : {
1677 28 : ret.push_back(array->GetName());
1678 : }
1679 3 : return ret;
1680 : }
1681 :
1682 : /************************************************************************/
1683 : /* OpenMDArray() */
1684 : /************************************************************************/
1685 :
1686 3 : std::shared_ptr<GDALMDArray> GRIBGroup::OpenMDArray(const std::string &osName,
1687 : CSLConstList) const
1688 : {
1689 12 : for (const auto &array : m_poArrays)
1690 : {
1691 11 : if (array->GetName() == osName)
1692 2 : return array;
1693 : }
1694 1 : return nullptr;
1695 : }
1696 :
1697 : /************************************************************************/
1698 : /* GRIBArray() */
1699 : /************************************************************************/
1700 :
1701 23 : GRIBArray::GRIBArray(const std::string &osName,
1702 23 : const std::shared_ptr<GRIBSharedResource> &poShared)
1703 : : GDALAbstractMDArray("/", osName),
1704 23 : GDALPamMDArray("/", osName, poShared->GetPAM()), m_poShared(poShared)
1705 : {
1706 23 : }
1707 :
1708 : /************************************************************************/
1709 : /* Init() */
1710 : /************************************************************************/
1711 :
1712 23 : void GRIBArray::Init(GRIBGroup *poGroup, GRIBDataset *poDS,
1713 : GRIBRasterBand *poBand, inventoryType *psInv)
1714 : {
1715 23 : std::shared_ptr<GDALDimension> poDimX;
1716 23 : std::shared_ptr<GDALDimension> poDimY;
1717 :
1718 : double adfGT[6];
1719 23 : poDS->GetGeoTransform(adfGT);
1720 :
1721 40 : for (int i = 1; i <= poGroup->m_nHorizDimCounter; i++)
1722 : {
1723 34 : std::string osXLookup("X");
1724 34 : std::string osYLookup("Y");
1725 34 : if (i > 1)
1726 : {
1727 15 : osXLookup += CPLSPrintf("%d", i);
1728 15 : osYLookup += CPLSPrintf("%d", i);
1729 : }
1730 34 : auto oIterX = poGroup->m_oMapDims.find(osXLookup);
1731 34 : auto oIterY = poGroup->m_oMapDims.find(osYLookup);
1732 34 : CPLAssert(oIterX != poGroup->m_oMapDims.end());
1733 34 : CPLAssert(oIterY != poGroup->m_oMapDims.end());
1734 34 : if (oIterX->second->GetSize() ==
1735 51 : static_cast<size_t>(poDS->GetRasterXSize()) &&
1736 17 : oIterY->second->GetSize() ==
1737 17 : static_cast<size_t>(poDS->GetRasterYSize()))
1738 : {
1739 17 : bool bOK = true;
1740 17 : auto poVar = oIterX->second->GetIndexingVariable();
1741 17 : constexpr double EPSILON = 1e-10;
1742 17 : if (poVar)
1743 : {
1744 17 : GUInt64 nStart = 0;
1745 17 : size_t nCount = 1;
1746 17 : double dfVal = 0;
1747 17 : poVar->Read(&nStart, &nCount, nullptr, nullptr, m_dt, &dfVal);
1748 17 : if (std::fabs(dfVal - (adfGT[0] + 0.5 * adfGT[1])) >
1749 17 : EPSILON * std::fabs(dfVal))
1750 : {
1751 0 : bOK = false;
1752 : }
1753 : }
1754 17 : if (bOK)
1755 : {
1756 17 : poVar = oIterY->second->GetIndexingVariable();
1757 17 : if (poVar)
1758 : {
1759 17 : GUInt64 nStart = 0;
1760 17 : size_t nCount = 1;
1761 17 : double dfVal = 0;
1762 17 : poVar->Read(&nStart, &nCount, nullptr, nullptr, m_dt,
1763 : &dfVal);
1764 17 : if (std::fabs(dfVal -
1765 17 : (adfGT[3] + poDS->nRasterYSize * adfGT[5] -
1766 17 : 0.5 * adfGT[5])) >
1767 17 : EPSILON * std::fabs(dfVal))
1768 : {
1769 0 : bOK = false;
1770 : }
1771 : }
1772 : }
1773 17 : if (bOK)
1774 : {
1775 17 : poDimX = oIterX->second;
1776 17 : poDimY = oIterY->second;
1777 17 : break;
1778 : }
1779 : }
1780 : }
1781 :
1782 23 : if (!poDimX || !poDimY)
1783 : {
1784 6 : poGroup->m_nHorizDimCounter++;
1785 : {
1786 12 : std::string osName("Y");
1787 6 : if (poGroup->m_nHorizDimCounter >= 2)
1788 : {
1789 2 : osName = CPLSPrintf("Y%d", poGroup->m_nHorizDimCounter);
1790 : }
1791 :
1792 12 : poDimY = std::make_shared<GDALDimensionWeakIndexingVar>(
1793 6 : poGroup->GetFullName(), osName, GDAL_DIM_TYPE_HORIZONTAL_Y,
1794 18 : std::string(), poDS->GetRasterYSize());
1795 6 : poGroup->m_oMapDims[osName] = poDimY;
1796 6 : poGroup->m_dims.emplace_back(poDimY);
1797 :
1798 : auto var = GDALMDArrayRegularlySpaced::Create(
1799 : "/", poDimY->GetName(), poDimY,
1800 12 : adfGT[3] + poDS->GetRasterYSize() * adfGT[5], -adfGT[5], 0.5);
1801 6 : poDimY->SetIndexingVariable(var);
1802 6 : poGroup->AddArray(var);
1803 : }
1804 :
1805 : {
1806 12 : std::string osName("X");
1807 6 : if (poGroup->m_nHorizDimCounter >= 2)
1808 : {
1809 2 : osName = CPLSPrintf("X%d", poGroup->m_nHorizDimCounter);
1810 : }
1811 :
1812 12 : poDimX = std::make_shared<GDALDimensionWeakIndexingVar>(
1813 6 : poGroup->GetFullName(), osName, GDAL_DIM_TYPE_HORIZONTAL_X,
1814 18 : std::string(), poDS->GetRasterXSize());
1815 6 : poGroup->m_oMapDims[osName] = poDimX;
1816 6 : poGroup->m_dims.emplace_back(poDimX);
1817 :
1818 : auto var = GDALMDArrayRegularlySpaced::Create(
1819 12 : "/", poDimX->GetName(), poDimX, adfGT[0], adfGT[1], 0.5);
1820 6 : poDimX->SetIndexingVariable(var);
1821 6 : poGroup->AddArray(var);
1822 : }
1823 : }
1824 :
1825 23 : m_dims.emplace_back(poDimY);
1826 23 : m_dims.emplace_back(poDimX);
1827 23 : if (poDS->m_poSRS.get())
1828 : {
1829 23 : m_poSRS.reset(poDS->m_poSRS->Clone());
1830 23 : if (poDS->m_poSRS->GetDataAxisToSRSAxisMapping() ==
1831 46 : std::vector<int>{2, 1})
1832 22 : m_poSRS->SetDataAxisToSRSAxisMapping({1, 2});
1833 : else
1834 1 : m_poSRS->SetDataAxisToSRSAxisMapping({2, 1});
1835 : }
1836 :
1837 : const char *pszGribNormalizeUnits =
1838 23 : CPLGetConfigOption("GRIB_NORMALIZE_UNITS", "YES");
1839 23 : const bool bMetricUnits = CPLTestBool(pszGribNormalizeUnits);
1840 :
1841 46 : m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
1842 46 : GetFullName(), "name", psInv->element));
1843 46 : m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
1844 23 : GetFullName(), "long_name",
1845 69 : ConvertUnitInText(bMetricUnits, psInv->comment)));
1846 23 : m_osUnit = ConvertUnitInText(bMetricUnits, psInv->unitName);
1847 23 : if (!m_osUnit.empty() && m_osUnit[0] == '[' && m_osUnit.back() == ']')
1848 : {
1849 23 : m_osUnit = m_osUnit.substr(1, m_osUnit.size() - 2);
1850 : }
1851 46 : m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
1852 46 : GetFullName(), "first_level", psInv->shortFstLevel));
1853 :
1854 23 : if (poBand->m_nDisciplineCode >= 0)
1855 : {
1856 14 : m_attributes.emplace_back(std::make_shared<GDALAttributeNumeric>(
1857 14 : GetFullName(), "discipline_code", poBand->m_nDisciplineCode));
1858 : }
1859 23 : if (!poBand->m_osDisciplineName.empty())
1860 : {
1861 14 : m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
1862 14 : GetFullName(), "discipline_name", poBand->m_osDisciplineName));
1863 : }
1864 23 : if (poBand->m_nCenter >= 0)
1865 : {
1866 14 : m_attributes.emplace_back(std::make_shared<GDALAttributeNumeric>(
1867 14 : GetFullName(), "center_code", poBand->m_nCenter));
1868 : }
1869 23 : if (!poBand->m_osCenterName.empty())
1870 : {
1871 14 : m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
1872 14 : GetFullName(), "center_name", poBand->m_osCenterName));
1873 : }
1874 23 : if (poBand->m_nSubCenter >= 0)
1875 : {
1876 12 : m_attributes.emplace_back(std::make_shared<GDALAttributeNumeric>(
1877 12 : GetFullName(), "subcenter_code", poBand->m_nSubCenter));
1878 : }
1879 23 : if (!poBand->m_osSubCenterName.empty())
1880 : {
1881 0 : m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
1882 0 : GetFullName(), "subcenter_name", poBand->m_osSubCenterName));
1883 : }
1884 23 : if (!poBand->m_osSignRefTimeName.empty())
1885 : {
1886 14 : m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
1887 7 : GetFullName(), "signification_of_ref_time",
1888 14 : poBand->m_osSignRefTimeName));
1889 : }
1890 23 : if (!poBand->m_osRefTime.empty())
1891 : {
1892 14 : m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
1893 14 : GetFullName(), "reference_time_iso8601", poBand->m_osRefTime));
1894 : }
1895 23 : if (!poBand->m_osProductionStatus.empty())
1896 : {
1897 14 : m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
1898 14 : GetFullName(), "production_status", poBand->m_osProductionStatus));
1899 : }
1900 23 : if (!poBand->m_osType.empty())
1901 : {
1902 14 : m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
1903 14 : GetFullName(), "type", poBand->m_osType));
1904 : }
1905 23 : if (poBand->m_nPDTN >= 0)
1906 : {
1907 14 : m_attributes.emplace_back(std::make_shared<GDALAttributeNumeric>(
1908 7 : GetFullName(), "product_definition_template_number",
1909 14 : poBand->m_nPDTN));
1910 : }
1911 23 : if (!poBand->m_anPDSTemplateAssembledValues.empty())
1912 : {
1913 14 : m_attributes.emplace_back(std::make_shared<GDALAttributeNumeric>(
1914 7 : GetFullName(), "product_definition_numbers",
1915 14 : poBand->m_anPDSTemplateAssembledValues));
1916 : }
1917 :
1918 23 : int bHasNoData = FALSE;
1919 23 : double dfNoData = poBand->GetNoDataValue(&bHasNoData);
1920 23 : if (bHasNoData)
1921 : {
1922 14 : m_abyNoData.resize(sizeof(double));
1923 14 : memcpy(&m_abyNoData[0], &dfNoData, sizeof(double));
1924 : }
1925 23 : }
1926 :
1927 : /************************************************************************/
1928 : /* ExtendTimeDim() */
1929 : /************************************************************************/
1930 :
1931 24 : void GRIBArray::ExtendTimeDim(vsi_l_offset nOffset, int subgNum,
1932 : double dfValidTime)
1933 : {
1934 24 : m_anOffsets.push_back(nOffset);
1935 24 : m_anSubgNums.push_back(subgNum);
1936 24 : m_adfTimes.push_back(dfValidTime);
1937 24 : }
1938 :
1939 : /************************************************************************/
1940 : /* Finalize() */
1941 : /************************************************************************/
1942 :
1943 23 : void GRIBArray::Finalize(GRIBGroup *poGroup, inventoryType *psInv)
1944 : {
1945 23 : CPLAssert(!m_adfTimes.empty());
1946 23 : CPLAssert(m_anOffsets.size() == m_adfTimes.size());
1947 :
1948 23 : if (m_adfTimes.size() == 1)
1949 : {
1950 44 : m_attributes.emplace_back(std::make_shared<GDALAttributeNumeric>(
1951 44 : GetFullName(), "forecast_time", psInv->foreSec));
1952 44 : m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
1953 44 : GetFullName(), "forecast_time_unit", "sec"));
1954 44 : m_attributes.emplace_back(std::make_shared<GDALAttributeNumeric>(
1955 44 : GetFullName(), "reference_time", psInv->refTime));
1956 44 : m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
1957 44 : GetFullName(), "reference_time_unit", "sec UTC"));
1958 44 : m_attributes.emplace_back(std::make_shared<GDALAttributeNumeric>(
1959 44 : GetFullName(), "validity_time", m_adfTimes[0]));
1960 44 : m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
1961 44 : GetFullName(), "validity_time_unit", "sec UTC"));
1962 22 : return;
1963 : }
1964 :
1965 1 : std::shared_ptr<GDALDimension> poDimTime;
1966 :
1967 3 : for (const auto &poDim : poGroup->m_dims)
1968 : {
1969 2 : if (STARTS_WITH(poDim->GetName().c_str(), "TIME") &&
1970 0 : poDim->GetSize() == m_adfTimes.size())
1971 : {
1972 0 : auto poVar = poDim->GetIndexingVariable();
1973 0 : if (poVar)
1974 : {
1975 0 : GUInt64 nStart = 0;
1976 0 : size_t nCount = 1;
1977 0 : double dfStartTime = 0;
1978 0 : poVar->Read(&nStart, &nCount, nullptr, nullptr, m_dt,
1979 : &dfStartTime);
1980 0 : if (dfStartTime == m_adfTimes[0])
1981 : {
1982 0 : poDimTime = poDim;
1983 0 : break;
1984 : }
1985 : }
1986 : }
1987 : }
1988 :
1989 1 : if (!poDimTime)
1990 : {
1991 2 : std::string osName("TIME");
1992 1 : int counter = 2;
1993 1 : while (poGroup->m_oMapDims.find(osName) != poGroup->m_oMapDims.end())
1994 : {
1995 0 : osName = CPLSPrintf("TIME%d", counter);
1996 0 : counter++;
1997 : }
1998 :
1999 2 : poDimTime = std::make_shared<GDALDimensionWeakIndexingVar>(
2000 1 : poGroup->GetFullName(), osName, GDAL_DIM_TYPE_TEMPORAL,
2001 3 : std::string(), m_adfTimes.size());
2002 1 : poGroup->m_oMapDims[osName] = poDimTime;
2003 1 : poGroup->m_dims.emplace_back(poDimTime);
2004 :
2005 1 : auto var = poGroup->m_memRootGroup->CreateMDArray(
2006 : poDimTime->GetName(),
2007 4 : std::vector<std::shared_ptr<GDALDimension>>{poDimTime},
2008 3 : GDALExtendedDataType::Create(GDT_Float64), nullptr);
2009 1 : poDimTime->SetIndexingVariable(var);
2010 1 : poGroup->AddArray(var);
2011 :
2012 1 : GUInt64 nStart = 0;
2013 1 : size_t nCount = m_adfTimes.size();
2014 1 : var->SetUnit("sec UTC");
2015 1 : const GUInt64 anStart[] = {nStart};
2016 1 : const size_t anCount[] = {nCount};
2017 2 : var->Write(anStart, anCount, nullptr, nullptr, var->GetDataType(),
2018 1 : &m_adfTimes[0]);
2019 1 : auto attr = var->CreateAttribute("long_name", {},
2020 3 : GDALExtendedDataType::CreateString());
2021 1 : attr->Write("validity_time");
2022 : }
2023 :
2024 1 : m_dims.insert(m_dims.begin(), poDimTime);
2025 1 : if (m_poSRS)
2026 : {
2027 2 : auto mapping = m_poSRS->GetDataAxisToSRSAxisMapping();
2028 3 : for (auto &v : mapping)
2029 2 : v += 1;
2030 1 : m_poSRS->SetDataAxisToSRSAxisMapping(mapping);
2031 : }
2032 : }
2033 :
2034 : /************************************************************************/
2035 : /* LoadData() */
2036 : /************************************************************************/
2037 :
2038 6 : const std::vector<double> &GRIBSharedResource::LoadData(vsi_l_offset nOffset,
2039 : int subgNum)
2040 : {
2041 6 : if (m_nOffsetCurData == nOffset)
2042 1 : return m_adfCurData;
2043 :
2044 5 : grib_MetaData *metadata = nullptr;
2045 5 : double *data = nullptr;
2046 5 : GRIBRasterBand::ReadGribData(m_fp, nOffset, subgNum, &data, &metadata);
2047 5 : if (data == nullptr || metadata == nullptr)
2048 : {
2049 0 : if (metadata != nullptr)
2050 : {
2051 0 : MetaFree(metadata);
2052 0 : delete metadata;
2053 : }
2054 0 : free(data);
2055 0 : m_adfCurData.clear();
2056 0 : return m_adfCurData;
2057 : }
2058 5 : const int nx = metadata->gds.Nx;
2059 5 : const int ny = metadata->gds.Ny;
2060 5 : MetaFree(metadata);
2061 5 : delete metadata;
2062 5 : if (nx <= 0 || ny <= 0)
2063 : {
2064 0 : free(data);
2065 0 : m_adfCurData.clear();
2066 0 : return m_adfCurData;
2067 : }
2068 5 : const size_t nPointCount = static_cast<size_t>(nx) * ny;
2069 5 : const size_t nByteCount = nPointCount * sizeof(double);
2070 : #ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
2071 : if (nByteCount > static_cast<size_t>(INT_MAX))
2072 : {
2073 : CPLError(CE_Failure, CPLE_OutOfMemory,
2074 : "Too large memory allocation attempt");
2075 : free(data);
2076 : m_adfCurData.clear();
2077 : return m_adfCurData;
2078 : }
2079 : #endif
2080 : try
2081 : {
2082 5 : m_adfCurData.resize(nPointCount);
2083 : }
2084 0 : catch (const std::exception &e)
2085 : {
2086 0 : CPLError(CE_Failure, CPLE_OutOfMemory, "%s", e.what());
2087 0 : free(data);
2088 0 : m_adfCurData.clear();
2089 0 : return m_adfCurData;
2090 : }
2091 5 : m_nOffsetCurData = nOffset;
2092 5 : memcpy(&m_adfCurData[0], data, nByteCount);
2093 5 : free(data);
2094 5 : return m_adfCurData;
2095 : }
2096 :
2097 : /************************************************************************/
2098 : /* IRead() */
2099 : /************************************************************************/
2100 :
2101 4 : bool GRIBArray::IRead(const GUInt64 *arrayStartIdx, const size_t *count,
2102 : const GInt64 *arrayStep, const GPtrDiff_t *bufferStride,
2103 : const GDALExtendedDataType &bufferDataType,
2104 : void *pDstBuffer) const
2105 : {
2106 4 : const size_t nBufferDTSize(bufferDataType.GetSize());
2107 4 : if (m_dims.size() == 2)
2108 : {
2109 2 : auto &vals = m_poShared->LoadData(m_anOffsets[0], m_anSubgNums[0]);
2110 2 : constexpr int Y_IDX = 0;
2111 2 : constexpr int X_IDX = 1;
2112 4 : if (vals.empty() ||
2113 2 : vals.size() != m_dims[Y_IDX]->GetSize() * m_dims[X_IDX]->GetSize())
2114 0 : return false;
2115 2 : const size_t nWidth(static_cast<size_t>(m_dims[X_IDX]->GetSize()));
2116 2 : const bool bDirectCopy = m_dt == bufferDataType &&
2117 3 : arrayStep[X_IDX] == 1 &&
2118 1 : bufferStride[X_IDX] == 1;
2119 150 : for (size_t j = 0; j < count[Y_IDX]; j++)
2120 : {
2121 148 : const size_t y = static_cast<size_t>(arrayStartIdx[Y_IDX] +
2122 148 : j * arrayStep[Y_IDX]);
2123 148 : GByte *pabyDstPtr = static_cast<GByte *>(pDstBuffer) +
2124 148 : j * bufferStride[Y_IDX] * nBufferDTSize;
2125 148 : const size_t x = static_cast<size_t>(arrayStartIdx[X_IDX]);
2126 148 : const double *srcPtr = &vals[y * nWidth + x];
2127 148 : if (bDirectCopy)
2128 : {
2129 74 : memcpy(pabyDstPtr, srcPtr, count[X_IDX] * sizeof(double));
2130 : }
2131 : else
2132 : {
2133 74 : const auto dstPtrInc = bufferStride[X_IDX] * nBufferDTSize;
2134 4958 : for (size_t i = 0; i < count[X_IDX]; i++)
2135 : {
2136 4884 : GDALExtendedDataType::CopyValue(srcPtr, m_dt, pabyDstPtr,
2137 : bufferDataType);
2138 4884 : srcPtr += static_cast<std::ptrdiff_t>(arrayStep[X_IDX]);
2139 4884 : pabyDstPtr += dstPtrInc;
2140 : }
2141 : }
2142 : }
2143 2 : return true;
2144 : }
2145 :
2146 2 : constexpr int T_IDX = 0;
2147 2 : constexpr int Y_IDX = 1;
2148 2 : constexpr int X_IDX = 2;
2149 2 : const size_t nWidth(static_cast<size_t>(m_dims[X_IDX]->GetSize()));
2150 3 : const bool bDirectCopy = m_dt == bufferDataType && arrayStep[X_IDX] == 1 &&
2151 1 : bufferStride[X_IDX] == 1;
2152 6 : for (size_t k = 0; k < count[T_IDX]; k++)
2153 : {
2154 4 : const size_t tIdx =
2155 4 : static_cast<size_t>(arrayStartIdx[T_IDX] + k * arrayStep[T_IDX]);
2156 4 : CPLAssert(tIdx < m_anOffsets.size());
2157 : auto &vals =
2158 4 : m_poShared->LoadData(m_anOffsets[tIdx], m_anSubgNums[tIdx]);
2159 8 : if (vals.empty() ||
2160 4 : vals.size() != m_dims[Y_IDX]->GetSize() * m_dims[X_IDX]->GetSize())
2161 0 : return false;
2162 520 : for (size_t j = 0; j < count[Y_IDX]; j++)
2163 : {
2164 516 : const size_t y = static_cast<size_t>(arrayStartIdx[Y_IDX] +
2165 516 : j * arrayStep[Y_IDX]);
2166 516 : GByte *pabyDstPtr =
2167 : static_cast<GByte *>(pDstBuffer) +
2168 516 : (k * bufferStride[T_IDX] + j * bufferStride[Y_IDX]) *
2169 : nBufferDTSize;
2170 516 : const size_t x = static_cast<size_t>(arrayStartIdx[X_IDX]);
2171 516 : const double *srcPtr = &vals[y * nWidth + x];
2172 516 : if (bDirectCopy)
2173 : {
2174 258 : memcpy(pabyDstPtr, srcPtr, count[X_IDX] * sizeof(double));
2175 : }
2176 : else
2177 : {
2178 258 : const auto dstPtrInc = bufferStride[X_IDX] * nBufferDTSize;
2179 45924 : for (size_t i = 0; i < count[X_IDX]; i++)
2180 : {
2181 45666 : GDALExtendedDataType::CopyValue(srcPtr, m_dt, pabyDstPtr,
2182 : bufferDataType);
2183 45666 : srcPtr += static_cast<std::ptrdiff_t>(arrayStep[X_IDX]);
2184 45666 : pabyDstPtr += dstPtrInc;
2185 : }
2186 : }
2187 : }
2188 : }
2189 :
2190 2 : return true;
2191 : }
2192 :
2193 : /************************************************************************/
2194 : /* OpenMultiDim() */
2195 : /************************************************************************/
2196 :
2197 4 : GDALDataset *GRIBDataset::OpenMultiDim(GDALOpenInfo *poOpenInfo)
2198 :
2199 : {
2200 : auto poShared = std::make_shared<GRIBSharedResource>(
2201 8 : poOpenInfo->pszFilename, poOpenInfo->fpL);
2202 8 : auto poRootGroup = std::make_shared<GRIBGroup>(poShared);
2203 4 : poOpenInfo->fpL = nullptr;
2204 :
2205 4 : VSIFSeekL(poShared->m_fp, 0, SEEK_SET);
2206 :
2207 : // Contains an GRIB2 message inventory of the file.
2208 : // We can't use the potential .idx file
2209 8 : auto pInventories = std::make_unique<InventoryWrapperGrib>(poShared->m_fp);
2210 :
2211 4 : if (pInventories->result() <= 0)
2212 : {
2213 0 : char *errMsg = errSprintf(nullptr);
2214 0 : if (errMsg != nullptr)
2215 0 : CPLDebug("GRIB", "%s", errMsg);
2216 0 : free(errMsg);
2217 :
2218 0 : CPLError(CE_Failure, CPLE_OpenFailed,
2219 : "%s is a grib file, "
2220 : "but no raster dataset was successfully identified.",
2221 : poOpenInfo->pszFilename);
2222 0 : return nullptr;
2223 : }
2224 :
2225 : // Create band objects.
2226 4 : GRIBDataset *poDS = new GRIBDataset();
2227 4 : poDS->fp = poShared->m_fp;
2228 :
2229 4 : std::shared_ptr<GRIBArray> poArray;
2230 8 : std::set<std::string> oSetArrayNames;
2231 8 : std::string osElement;
2232 8 : std::string osShortFstLevel;
2233 4 : double dfRefTime = 0;
2234 4 : double dfForecastTime = 0;
2235 28 : for (uInt4 i = 0; i < pInventories->length(); ++i)
2236 : {
2237 24 : inventoryType *psInv = pInventories->get(i);
2238 24 : uInt4 bandNr = i + 1;
2239 :
2240 : // GRIB messages can be preceded by "garbage". GRIB2Inventory()
2241 : // does not return the offset to the real start of the message
2242 : GByte abyHeader[1024 + 1];
2243 24 : VSIFSeekL(poShared->m_fp, psInv->start, SEEK_SET);
2244 : size_t nRead =
2245 24 : VSIFReadL(abyHeader, 1, sizeof(abyHeader) - 1, poShared->m_fp);
2246 24 : abyHeader[nRead] = 0;
2247 : // Find the real offset of the fist message
2248 24 : const char *pasHeader = reinterpret_cast<char *>(abyHeader);
2249 24 : int nOffsetFirstMessage = 0;
2250 144 : for (int j = 0; j < poOpenInfo->nHeaderBytes - 3; j++)
2251 : {
2252 144 : if (STARTS_WITH_CI(pasHeader + j, "GRIB")
2253 : #ifdef ENABLE_TDLP
2254 : || STARTS_WITH_CI(pasHeader + j, "TDLP")
2255 : #endif
2256 : )
2257 : {
2258 24 : nOffsetFirstMessage = j;
2259 24 : break;
2260 : }
2261 : }
2262 24 : psInv->start += nOffsetFirstMessage;
2263 :
2264 44 : if (poArray == nullptr || osElement != psInv->element ||
2265 45 : osShortFstLevel != psInv->shortFstLevel ||
2266 1 : !((dfRefTime == psInv->refTime &&
2267 1 : dfForecastTime != psInv->foreSec) ||
2268 0 : (dfRefTime != psInv->refTime &&
2269 0 : dfForecastTime == psInv->foreSec)))
2270 : {
2271 23 : if (poArray)
2272 : {
2273 19 : poArray->Finalize(poRootGroup.get(), pInventories->get(i - 1));
2274 19 : poRootGroup->AddArray(poArray);
2275 : }
2276 :
2277 23 : grib_MetaData *metaData = nullptr;
2278 23 : GRIBRasterBand::ReadGribData(poShared->m_fp, psInv->start,
2279 23 : psInv->subgNum, nullptr, &metaData);
2280 23 : if (metaData == nullptr || metaData->gds.Nx < 1 ||
2281 23 : metaData->gds.Ny < 1)
2282 : {
2283 0 : CPLError(CE_Failure, CPLE_OpenFailed,
2284 : "%s is a grib file, "
2285 : "but no raster dataset was successfully identified.",
2286 : poOpenInfo->pszFilename);
2287 :
2288 0 : if (metaData != nullptr)
2289 : {
2290 0 : MetaFree(metaData);
2291 0 : delete metaData;
2292 : }
2293 0 : poDS->fp = nullptr;
2294 0 : CPLReleaseMutex(hGRIBMutex);
2295 0 : delete poDS;
2296 0 : CPLAcquireMutex(hGRIBMutex, 1000.0);
2297 0 : return nullptr;
2298 : }
2299 23 : psInv->GribVersion = metaData->GribVersion;
2300 :
2301 : // Set the DataSet's x,y size, georeference and projection from
2302 : // the first GRIB band.
2303 23 : poDS->SetGribMetaData(metaData);
2304 :
2305 : // coverity[tainted_data]
2306 46 : GRIBRasterBand gribBand(poDS, bandNr, psInv);
2307 23 : if (psInv->GribVersion == 2)
2308 7 : gribBand.FindPDSTemplateGRIB2();
2309 23 : osElement = psInv->element;
2310 23 : osShortFstLevel = psInv->shortFstLevel;
2311 23 : dfRefTime = psInv->refTime;
2312 23 : dfForecastTime = psInv->foreSec;
2313 :
2314 46 : std::string osRadix(osElement + '_' + osShortFstLevel);
2315 46 : std::string osName(osRadix);
2316 23 : int nCounter = 2;
2317 23 : while (oSetArrayNames.find(osName) != oSetArrayNames.end())
2318 : {
2319 0 : osName = osRadix + CPLSPrintf("_%d", nCounter);
2320 0 : nCounter++;
2321 : }
2322 23 : oSetArrayNames.insert(osName);
2323 23 : poArray = GRIBArray::Create(osName, poShared);
2324 23 : poArray->Init(poRootGroup.get(), poDS, &gribBand, psInv);
2325 23 : poArray->ExtendTimeDim(psInv->start, psInv->subgNum,
2326 : psInv->validTime);
2327 :
2328 23 : MetaFree(metaData);
2329 23 : delete metaData;
2330 : }
2331 : else
2332 : {
2333 1 : poArray->ExtendTimeDim(psInv->start, psInv->subgNum,
2334 : psInv->validTime);
2335 : }
2336 : }
2337 :
2338 4 : if (poArray)
2339 : {
2340 8 : poArray->Finalize(poRootGroup.get(),
2341 4 : pInventories->get(pInventories->length() - 1));
2342 4 : poRootGroup->AddArray(poArray);
2343 : }
2344 :
2345 4 : poDS->fp = nullptr;
2346 4 : poDS->m_poRootGroup = poRootGroup;
2347 :
2348 4 : poDS->SetDescription(poOpenInfo->pszFilename);
2349 :
2350 : // Release hGRIBMutex otherwise we'll deadlock with GDALDataset own
2351 : // hGRIBMutex.
2352 4 : CPLReleaseMutex(hGRIBMutex);
2353 4 : poDS->TryLoadXML();
2354 4 : CPLAcquireMutex(hGRIBMutex, 1000.0);
2355 :
2356 4 : return poDS;
2357 : }
2358 :
2359 : /************************************************************************/
2360 : /* SetMetadata() */
2361 : /************************************************************************/
2362 :
2363 317 : void GRIBDataset::SetGribMetaData(grib_MetaData *meta)
2364 : {
2365 317 : nRasterXSize = meta->gds.Nx;
2366 317 : nRasterYSize = meta->gds.Ny;
2367 :
2368 : // Image projection.
2369 317 : OGRSpatialReference oSRS;
2370 317 : oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2371 :
2372 317 : switch (meta->gds.projType)
2373 : {
2374 128 : case GS3_LATLON:
2375 : case GS3_GAUSSIAN_LATLON:
2376 : // No projection, only latlon system (geographic).
2377 128 : break;
2378 2 : case GS3_ROTATED_LATLON:
2379 : // will apply pole rotation afterwards
2380 2 : break;
2381 26 : case GS3_MERCATOR:
2382 26 : if (meta->gds.orientLon == 0.0)
2383 : {
2384 26 : if (meta->gds.meshLat == 0.0)
2385 5 : oSRS.SetMercator(0.0, 0.0, 1.0, 0.0, 0.0);
2386 : else
2387 21 : oSRS.SetMercator2SP(meta->gds.meshLat, 0.0, 0.0, 0.0, 0.0);
2388 : }
2389 : else
2390 : {
2391 0 : CPLError(CE_Warning, CPLE_NotSupported,
2392 : "Orientation of the grid != 0 not supported");
2393 0 : return;
2394 : }
2395 26 : break;
2396 133 : case GS3_TRANSVERSE_MERCATOR:
2397 133 : oSRS.SetTM(meta->gds.latitude_of_origin,
2398 : Lon360to180(meta->gds.central_meridian),
2399 133 : std::abs(meta->gds.scaleLat1 - 0.9996) < 1e8
2400 : ? 0.9996
2401 : : meta->gds.scaleLat1,
2402 : meta->gds.x0, meta->gds.y0);
2403 133 : break;
2404 9 : case GS3_POLAR:
2405 9 : oSRS.SetPS(meta->gds.meshLat, meta->gds.orientLon, 1.0, 0.0, 0.0);
2406 9 : break;
2407 8 : case GS3_LAMBERT:
2408 8 : oSRS.SetLCC(meta->gds.scaleLat1, meta->gds.scaleLat2,
2409 : meta->gds.meshLat, Lon360to180(meta->gds.orientLon),
2410 : 0.0, 0.0);
2411 8 : break;
2412 5 : case GS3_ALBERS_EQUAL_AREA:
2413 5 : oSRS.SetACEA(meta->gds.scaleLat1, meta->gds.scaleLat2,
2414 : meta->gds.meshLat, Lon360to180(meta->gds.orientLon),
2415 : 0.0, 0.0);
2416 5 : break;
2417 :
2418 0 : case GS3_ORTHOGRAPHIC:
2419 :
2420 : // oSRS.SetOrthographic( 0.0, meta->gds.orientLon,
2421 : // meta->gds.lon2, meta->gds.lat2);
2422 :
2423 : // oSRS.SetGEOS( meta->gds.orientLon, meta->gds.stretchFactor,
2424 : // meta->gds.lon2, meta->gds.lat2);
2425 :
2426 : // TODO: Hardcoded for now. How to parse the meta->gds section?
2427 0 : oSRS.SetGEOS(0, 35785831, 0, 0);
2428 0 : break;
2429 6 : case GS3_LAMBERT_AZIMUTHAL:
2430 6 : oSRS.SetLAEA(meta->gds.meshLat, Lon360to180(meta->gds.orientLon),
2431 : 0.0, 0.0);
2432 6 : break;
2433 :
2434 0 : case GS3_EQUATOR_EQUIDIST:
2435 0 : break;
2436 0 : case GS3_AZIMUTH_RANGE:
2437 0 : break;
2438 : }
2439 :
2440 317 : if (oSRS.IsProjected())
2441 : {
2442 187 : oSRS.SetLinearUnits("Metre", 1.0);
2443 : }
2444 :
2445 317 : const bool bHaveEarthModel =
2446 317 : meta->gds.majEarth > 0.0 && meta->gds.minEarth > 0.0;
2447 : // In meters.
2448 : const double a = bHaveEarthModel
2449 317 : ? meta->gds.majEarth * 1.0e3
2450 2 : : CPLAtof(CPLGetConfigOption("GRIB_DEFAULT_SEMI_MAJOR",
2451 317 : "6377563.396"));
2452 : const double b =
2453 : bHaveEarthModel
2454 319 : ? meta->gds.minEarth * 1.0e3
2455 2 : : (meta->gds.f_sphere
2456 2 : ? a
2457 0 : : CPLAtof(CPLGetConfigOption("GRIB_DEFAULT_SEMI_MINOR",
2458 317 : "6356256.910")));
2459 317 : if (meta->gds.majEarth == 0 || meta->gds.minEarth == 0)
2460 : {
2461 0 : CPLDebug("GRIB", "No earth model. Assuming a=%f and b=%f", a, b);
2462 : }
2463 317 : else if (meta->gds.majEarth < 0 || meta->gds.minEarth < 0)
2464 : {
2465 : const char *pszUseDefaultSpheroid =
2466 2 : CPLGetConfigOption("GRIB_USE_DEFAULT_SPHEROID", nullptr);
2467 2 : if (!pszUseDefaultSpheroid)
2468 : {
2469 1 : CPLError(CE_Failure, CPLE_AppDefined,
2470 : "The GRIB file contains invalid values for the spheroid. "
2471 : "You may set the GRIB_USE_DEFAULT_SPHEROID configuration "
2472 : "option to YES to use a default spheroid with "
2473 : "a=%f and b=%f",
2474 : a, b);
2475 1 : return;
2476 : }
2477 1 : else if (!CPLTestBool(pszUseDefaultSpheroid))
2478 : {
2479 0 : return;
2480 : }
2481 1 : CPLDebug("GRIB", "Invalid earth model. Assuming a=%f and b=%f", a, b);
2482 : }
2483 :
2484 316 : if (meta->gds.f_sphere || (a == b))
2485 : {
2486 81 : oSRS.SetGeogCS("Coordinate System imported from GRIB file", nullptr,
2487 : "Sphere", a, 0.0);
2488 : }
2489 : else
2490 : {
2491 235 : const double fInv = a / (a - b);
2492 335 : if (std::abs(a - 6378137.0) < 0.01 &&
2493 100 : std::abs(fInv - 298.257223563) < 1e-9) // WGS84
2494 : {
2495 98 : if (meta->gds.projType == GS3_LATLON)
2496 68 : oSRS.SetFromUserInput(SRS_WKT_WGS84_LAT_LONG);
2497 : else
2498 : {
2499 30 : oSRS.SetGeogCS("Coordinate System imported from GRIB file",
2500 : "WGS_1984", "WGS 84", 6378137., 298.257223563);
2501 : }
2502 : }
2503 139 : else if (std::abs(a - 6378137.0) < 0.01 &&
2504 2 : std::abs(fInv - 298.257222101) < 1e-9) // GRS80
2505 : {
2506 2 : oSRS.SetGeogCS("Coordinate System imported from GRIB file", nullptr,
2507 : "GRS80", 6378137., 298.257222101);
2508 : }
2509 : else
2510 : {
2511 135 : oSRS.SetGeogCS("Coordinate System imported from GRIB file", nullptr,
2512 : "Spheroid imported from GRIB file", a, fInv);
2513 : }
2514 : }
2515 :
2516 316 : if (meta->gds.projType == GS3_ROTATED_LATLON)
2517 : {
2518 2 : oSRS.SetDerivedGeogCRSWithPoleRotationGRIBConvention(
2519 : oSRS.GetName(), meta->gds.southLat, Lon360to180(meta->gds.southLon),
2520 : meta->gds.angleRotate);
2521 : }
2522 :
2523 316 : OGRSpatialReference oLL; // Construct the "geographic" part of oSRS.
2524 316 : oLL.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2525 316 : oLL.CopyGeogCSFrom(&oSRS);
2526 :
2527 316 : double rMinX = 0.0;
2528 316 : double rMaxY = 0.0;
2529 316 : double rPixelSizeX = 0.0;
2530 316 : double rPixelSizeY = 0.0;
2531 316 : bool bError = false;
2532 316 : if (meta->gds.projType == GS3_ORTHOGRAPHIC)
2533 : {
2534 : // This is what should work, but it doesn't. Dx seems to have an
2535 : // inverse relation with pixel size.
2536 : // rMinX = -meta->gds.Dx * (meta->gds.Nx / 2);
2537 : // rMaxY = meta->gds.Dy * (meta->gds.Ny / 2);
2538 : // Hardcoded for now, assumption: GEOS projection, full disc (like MSG).
2539 0 : const double geosExtentInMeters = 11137496.552;
2540 0 : rMinX = -(geosExtentInMeters / 2);
2541 0 : rMaxY = geosExtentInMeters / 2;
2542 0 : rPixelSizeX = geosExtentInMeters / meta->gds.Nx;
2543 0 : rPixelSizeY = geosExtentInMeters / meta->gds.Ny;
2544 : }
2545 316 : else if (meta->gds.projType == GS3_TRANSVERSE_MERCATOR)
2546 : {
2547 132 : rMinX = meta->gds.x1;
2548 132 : rMaxY = meta->gds.y2;
2549 132 : rPixelSizeX = meta->gds.Dx;
2550 132 : rPixelSizeY = meta->gds.Dy;
2551 : }
2552 184 : else if (oSRS.IsProjected() && meta->gds.projType != GS3_ROTATED_LATLON)
2553 : {
2554 : // Longitude in degrees, to be transformed to meters (or degrees in
2555 : // case of latlon).
2556 54 : rMinX = meta->gds.lon1;
2557 : // Latitude in degrees, to be transformed to meters.
2558 54 : rMaxY = meta->gds.lat1;
2559 :
2560 54 : if (m_poSRS == nullptr || m_poLL == nullptr ||
2561 54 : !m_poSRS->IsSame(&oSRS) || !m_poLL->IsSame(&oLL))
2562 : {
2563 108 : m_poCT = std::unique_ptr<OGRCoordinateTransformation>(
2564 54 : OGRCreateCoordinateTransformation(&(oLL), &(oSRS)));
2565 : }
2566 :
2567 : // Transform it to meters.
2568 54 : if ((m_poCT != nullptr) && m_poCT->Transform(1, &rMinX, &rMaxY))
2569 : {
2570 54 : if (meta->gds.scan == GRIB2BIT_2) // Y is minY, GDAL wants maxY.
2571 : {
2572 : // -1 because we GDAL needs the coordinates of the centre of
2573 : // the pixel.
2574 54 : rMaxY += (meta->gds.Ny - 1) * meta->gds.Dy;
2575 : }
2576 54 : rPixelSizeX = meta->gds.Dx;
2577 54 : rPixelSizeY = meta->gds.Dy;
2578 : }
2579 : else
2580 : {
2581 0 : rMinX = 0.0;
2582 0 : rMaxY = 0.0;
2583 :
2584 0 : rPixelSizeX = 1.0;
2585 0 : rPixelSizeY = -1.0;
2586 :
2587 0 : bError = true;
2588 0 : CPLError(CE_Warning, CPLE_AppDefined,
2589 : "Unable to perform coordinate transformations, so the "
2590 : "correct projected geotransform could not be deduced "
2591 : "from the lat/long control points. "
2592 : "Defaulting to ungeoreferenced.");
2593 : }
2594 : }
2595 : else
2596 : {
2597 : // Longitude in degrees, to be transformed to meters (or degrees in
2598 : // case of latlon).
2599 130 : rMinX = meta->gds.lon1;
2600 : // Latitude in degrees, to be transformed to meters.
2601 130 : rMaxY = meta->gds.lat1;
2602 :
2603 130 : double rMinY = meta->gds.lat2;
2604 130 : double rMaxX = meta->gds.lon2;
2605 130 : if (meta->gds.lat2 > rMaxY)
2606 : {
2607 100 : rMaxY = meta->gds.lat2;
2608 100 : rMinY = meta->gds.lat1;
2609 : }
2610 :
2611 130 : if (meta->gds.Nx == 1)
2612 16 : rPixelSizeX = meta->gds.Dx;
2613 114 : else if (meta->gds.lon1 > meta->gds.lon2)
2614 13 : rPixelSizeX = (360.0 - (meta->gds.lon1 - meta->gds.lon2)) /
2615 13 : (meta->gds.Nx - 1);
2616 : else
2617 101 : rPixelSizeX =
2618 101 : (meta->gds.lon2 - meta->gds.lon1) / (meta->gds.Nx - 1);
2619 :
2620 130 : if (meta->gds.Ny == 1)
2621 16 : rPixelSizeY = meta->gds.Dy;
2622 : else
2623 114 : rPixelSizeY = (rMaxY - rMinY) / (meta->gds.Ny - 1);
2624 :
2625 : // Do some sanity checks for cases that can't be handled by the above
2626 : // pixel size corrections. GRIB1 has a minimum precision of 0.001
2627 : // for latitudes and longitudes, so we'll allow a bit higher than that.
2628 130 : if (rPixelSizeX < 0 || fabs(rPixelSizeX - meta->gds.Dx) > 0.002)
2629 0 : rPixelSizeX = meta->gds.Dx;
2630 :
2631 130 : if (rPixelSizeY < 0 || fabs(rPixelSizeY - meta->gds.Dy) > 0.002)
2632 0 : rPixelSizeY = meta->gds.Dy;
2633 :
2634 : // GRIB2 files have longitudes in the [0-360] range
2635 : // Shift them to the traditional [-180,180] longitude range
2636 : // See https://github.com/OSGeo/gdal/issues/4524
2637 182 : if ((rMinX + rPixelSizeX >= 180 || rMaxX - rPixelSizeX >= 180) &&
2638 52 : CPLTestBool(
2639 : CPLGetConfigOption("GRIB_ADJUST_LONGITUDE_RANGE", "YES")))
2640 : {
2641 50 : if (rPixelSizeX * nRasterXSize > 360 + rPixelSizeX / 4)
2642 0 : CPLDebug("GRIB", "Cannot properly handle GRIB2 files with "
2643 : "overlaps and 0-360 longitudes");
2644 50 : else if (fabs(360 - rPixelSizeX * nRasterXSize) < rPixelSizeX / 4 &&
2645 19 : rMinX <= 180 && meta->gds.projType == GS3_LATLON)
2646 : {
2647 : // Find the first row number east of the antimeridian
2648 9 : const int nSplitAndSwapColumnCandidate =
2649 9 : static_cast<int>(ceil((180 - rMinX) / rPixelSizeX));
2650 9 : if (nSplitAndSwapColumnCandidate < nRasterXSize)
2651 : {
2652 9 : nSplitAndSwapColumn = nSplitAndSwapColumnCandidate;
2653 9 : CPLDebug("GRIB",
2654 : "Rewrapping around the antimeridian at column %d",
2655 : nSplitAndSwapColumn);
2656 9 : rMinX = -180;
2657 9 : }
2658 : }
2659 41 : else if (Lon360to180(rMinX) > Lon360to180(rMaxX))
2660 : {
2661 2 : CPLDebug("GRIB", "GRIB with 0-360 longitudes spanning across "
2662 : "the antimeridian");
2663 2 : rMinX = Lon360to180(rMinX);
2664 : }
2665 : else
2666 : {
2667 39 : CPLDebug("GRIB", "Shifting longitudes from %lf:%lf to %lf:%lf",
2668 : rMinX, rMaxX, Lon360to180(rMinX), Lon360to180(rMaxX));
2669 39 : rMinX = Lon360to180(rMinX);
2670 : }
2671 : }
2672 : }
2673 :
2674 : // https://gdal.org/user/raster_data_model.html :
2675 : // we need the top left corner of the top left pixel.
2676 : // At the moment we have the center of the pixel.
2677 316 : rMinX -= rPixelSizeX / 2;
2678 316 : rMaxY += rPixelSizeY / 2;
2679 :
2680 316 : adfGeoTransform[0] = rMinX;
2681 316 : adfGeoTransform[3] = rMaxY;
2682 316 : adfGeoTransform[1] = rPixelSizeX;
2683 316 : adfGeoTransform[5] = -rPixelSizeY;
2684 :
2685 316 : if (bError)
2686 0 : m_poSRS.reset();
2687 : else
2688 316 : m_poSRS.reset(oSRS.Clone());
2689 316 : m_poLL = std::unique_ptr<OGRSpatialReference>(oLL.Clone());
2690 : }
2691 :
2692 : /************************************************************************/
2693 : /* GDALDeregister_GRIB() */
2694 : /************************************************************************/
2695 :
2696 852 : static void GDALDeregister_GRIB(GDALDriver *)
2697 : {
2698 852 : if (hGRIBMutex != nullptr)
2699 : {
2700 1 : CPLDestroyMutex(hGRIBMutex);
2701 1 : hGRIBMutex = nullptr;
2702 : }
2703 852 : }
2704 :
2705 : /************************************************************************/
2706 : /* GDALGRIBDriver */
2707 : /************************************************************************/
2708 :
2709 : class GDALGRIBDriver : public GDALDriver
2710 : {
2711 : bool m_bHasFullInitMetadata = false;
2712 :
2713 : public:
2714 1217 : GDALGRIBDriver() = default;
2715 :
2716 : char **GetMetadata(const char *pszDomain = "") override;
2717 : const char *GetMetadataItem(const char *pszName,
2718 : const char *pszDomain) override;
2719 : };
2720 :
2721 : /************************************************************************/
2722 : /* GetMetadata() */
2723 : /************************************************************************/
2724 :
2725 272 : char **GDALGRIBDriver::GetMetadata(const char *pszDomain)
2726 : {
2727 272 : if (pszDomain == nullptr || EQUAL(pszDomain, ""))
2728 : {
2729 : // Defer until necessary the setting of the CreationOptionList
2730 : // to let a chance to JPEG2000 drivers to have been loaded.
2731 272 : if (!m_bHasFullInitMetadata)
2732 : {
2733 7 : m_bHasFullInitMetadata = true;
2734 :
2735 14 : std::vector<CPLString> aosJ2KDrivers;
2736 35 : for (size_t i = 0; i < CPL_ARRAYSIZE(apszJ2KDrivers); i++)
2737 : {
2738 28 : if (GDALGetDriverByName(apszJ2KDrivers[i]) != nullptr)
2739 : {
2740 14 : aosJ2KDrivers.push_back(apszJ2KDrivers[i]);
2741 : }
2742 : }
2743 :
2744 : CPLString osCreationOptionList(
2745 : "<CreationOptionList>"
2746 : " <Option name='DATA_ENCODING' type='string-select' "
2747 : "default='AUTO' "
2748 : "description='How data is encoded internally'>"
2749 : " <Value>AUTO</Value>"
2750 : " <Value>SIMPLE_PACKING</Value>"
2751 : " <Value>COMPLEX_PACKING</Value>"
2752 14 : " <Value>IEEE_FLOATING_POINT</Value>");
2753 7 : if (GDALGetDriverByName("PNG") != nullptr)
2754 7 : osCreationOptionList += " <Value>PNG</Value>";
2755 7 : if (!aosJ2KDrivers.empty())
2756 7 : osCreationOptionList += " <Value>JPEG2000</Value>";
2757 : osCreationOptionList +=
2758 : " </Option>"
2759 : " <Option name='NBITS' type='int' default='0' "
2760 : "description='Number of bits per value'/>"
2761 : " <Option name='DECIMAL_SCALE_FACTOR' type='int' default='0' "
2762 : "description='Value such that raw values are multiplied by "
2763 : "10^DECIMAL_SCALE_FACTOR before integer encoding'/>"
2764 : " <Option name='SPATIAL_DIFFERENCING_ORDER' type='int' "
2765 : "default='0' "
2766 7 : "description='Order of spatial difference' min='0' max='2'/>";
2767 7 : if (!aosJ2KDrivers.empty())
2768 : {
2769 : osCreationOptionList +=
2770 : " <Option name='COMPRESSION_RATIO' type='int' "
2771 : "default='1' min='1' max='100' "
2772 : "description='N:1 target compression ratio for JPEG2000'/>"
2773 : " <Option name='JPEG2000_DRIVER' type='string-select' "
2774 7 : "description='Explicitly select a JPEG2000 driver'>";
2775 21 : for (size_t i = 0; i < aosJ2KDrivers.size(); i++)
2776 : {
2777 : osCreationOptionList +=
2778 14 : " <Value>" + aosJ2KDrivers[i] + "</Value>";
2779 : }
2780 7 : osCreationOptionList += " </Option>";
2781 : }
2782 : osCreationOptionList +=
2783 : " <Option name='DISCIPLINE' type='int' "
2784 : "description='Discipline of the processed data'/>"
2785 : " <Option name='IDS' type='string' "
2786 : "description='String equivalent to the GRIB_IDS metadata "
2787 : "item'/>"
2788 : " <Option name='IDS_CENTER' type='int' "
2789 : "description='Originating/generating center'/>"
2790 : " <Option name='IDS_SUBCENTER' type='int' "
2791 : "description='Originating/generating subcenter'/>"
2792 : " <Option name='IDS_MASTER_TABLE' type='int' "
2793 : "description='GRIB master tables version number'/>"
2794 : " <Option name='IDS_SIGNF_REF_TIME' type='int' "
2795 : "description='Significance of Reference Time'/>"
2796 : " <Option name='IDS_REF_TIME' type='string' "
2797 : "description='Reference time as YYYY-MM-DDTHH:MM:SSZ'/>"
2798 : " <Option name='IDS_PROD_STATUS' type='int' "
2799 : "description='Production Status of Processed data'/>"
2800 : " <Option name='IDS_TYPE' type='int' "
2801 : "description='Type of processed data'/>"
2802 : " <Option name='PDS_PDTN' type='int' "
2803 : "description='Product Definition Template Number'/>"
2804 : " <Option name='PDS_TEMPLATE_NUMBERS' type='string' "
2805 : "description='Product definition template raw numbers'/>"
2806 : " <Option name='PDS_TEMPLATE_ASSEMBLED_VALUES' type='string' "
2807 : "description='Product definition template assembled values'/>"
2808 : " <Option name='INPUT_UNIT' type='string' "
2809 : "description='Unit of input values. Only for temperatures. C "
2810 : "or K'/>"
2811 : " <Option name='BAND_*' type='string' "
2812 : "description='Override options at band level'/>"
2813 7 : "</CreationOptionList>";
2814 :
2815 7 : SetMetadataItem(GDAL_DMD_CREATIONOPTIONLIST, osCreationOptionList);
2816 : }
2817 : }
2818 272 : return GDALDriver::GetMetadata(pszDomain);
2819 : }
2820 :
2821 : /************************************************************************/
2822 : /* GetMetadataItem() */
2823 : /************************************************************************/
2824 :
2825 51090 : const char *GDALGRIBDriver::GetMetadataItem(const char *pszName,
2826 : const char *pszDomain)
2827 : {
2828 51090 : if (pszDomain == nullptr || EQUAL(pszDomain, ""))
2829 : {
2830 : // Defer until necessary the setting of the CreationOptionList
2831 : // to let a chance to JPEG2000 drivers to have been loaded.
2832 51090 : if (EQUAL(pszName, GDAL_DMD_CREATIONOPTIONLIST))
2833 159 : GetMetadata();
2834 : }
2835 51090 : return GDALDriver::GetMetadataItem(pszName, pszDomain);
2836 : }
2837 :
2838 : /************************************************************************/
2839 : /* GDALRegister_GRIB() */
2840 : /************************************************************************/
2841 :
2842 1512 : void GDALRegister_GRIB()
2843 :
2844 : {
2845 1512 : if (GDALGetDriverByName(DRIVER_NAME) != nullptr)
2846 295 : return;
2847 :
2848 1217 : GDALDriver *poDriver = new GDALGRIBDriver();
2849 1217 : GRIBDriverSetCommonMetadata(poDriver);
2850 :
2851 1217 : poDriver->pfnOpen = GRIBDataset::Open;
2852 1217 : poDriver->pfnCreateCopy = GRIBDataset::CreateCopy;
2853 1217 : poDriver->pfnUnloadDriver = GDALDeregister_GRIB;
2854 :
2855 : #ifdef USE_AEC
2856 1217 : poDriver->SetMetadataItem("HAVE_AEC", "YES");
2857 : #endif
2858 :
2859 1217 : GetGDALDriverManager()->RegisterDriver(poDriver);
2860 : }
|
