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