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