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