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