1 /******************************************************************************
2 * $Id: gribdataset.cpp 28459 2015-02-12 13:48:21Z rouault $
3 *
4 * Project: GRIB Driver
5 * Purpose: GDALDataset driver for GRIB translator for read support
6 * Author: Bas Retsios, retsios@itc.nl
7 *
8 ******************************************************************************
9 * Copyright (c) 2007, ITC
10 * Copyright (c) 2008-2013, Even Rouault <even dot rouault at mines-paris dot org>
11 *
12 * Permission is hereby granted, free of charge, to any person obtaining a
13 * copy of this software and associated documentation files (the "Software"),
14 * to deal in the Software without restriction, including without limitation
15 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 * and/or sell copies of the Software, and to permit persons to whom the
17 * Software is furnished to do so, subject to the following conditions:
18 *
19 * The above copyright notice and this permission notice shall be included
20 * in all copies or substantial portions of the Software.
21 *
22 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 * DEALINGS IN THE SOFTWARE.
29 ******************************************************************************
30 *
31 */
32
33 #include "gdal_pam.h"
34 #include "cpl_multiproc.h"
35
36 #include "degrib18/degrib/degrib2.h"
37 #include "degrib18/degrib/inventory.h"
38 #include "degrib18/degrib/myerror.h"
39 #include "degrib18/degrib/filedatasource.h"
40 #include "degrib18/degrib/memorydatasource.h"
41
42 #include "ogr_spatialref.h"
43
44 CPL_CVSID("$Id: gribdataset.cpp 28459 2015-02-12 13:48:21Z rouault $");
45
46 CPL_C_START
47 void GDALRegister_GRIB(void);
48 CPL_C_END
49
50 static CPLMutex *hGRIBMutex = NULL;
51
52 /************************************************************************/
53 /* ==================================================================== */
54 /* GRIBDataset */
55 /* ==================================================================== */
56 /************************************************************************/
57
58 class GRIBRasterBand;
59
60 class GRIBDataset : public GDALPamDataset
61 {
62 friend class GRIBRasterBand;
63
64 public:
65 GRIBDataset();
66 ~GRIBDataset();
67
68 static GDALDataset *Open( GDALOpenInfo * );
69 static int Identify( GDALOpenInfo * );
70
71 CPLErr GetGeoTransform( double * padfTransform );
72 const char *GetProjectionRef();
73
74 private:
75 void SetGribMetaData(grib_MetaData* meta);
76 VSILFILE *fp;
77 char *pszProjection;
78 OGRCoordinateTransformation *poTransform;
79 double adfGeoTransform[6]; // Calculate and store once as GetGeoTransform may be called multiple times
80
81 GIntBig nCachedBytes;
82 GIntBig nCachedBytesThreshold;
83 int bCacheOnlyOneBand;
84 GRIBRasterBand* poLastUsedBand;
85 };
86
87 /************************************************************************/
88 /* ==================================================================== */
89 /* GRIBRasterBand */
90 /* ==================================================================== */
91 /************************************************************************/
92
93 class GRIBRasterBand : public GDALPamRasterBand
94 {
95 friend class GRIBDataset;
96
97 public:
98 GRIBRasterBand( GRIBDataset*, int, inventoryType* );
99 virtual ~GRIBRasterBand();
100 virtual CPLErr IReadBlock( int, int, void * );
101 virtual const char *GetDescription() const;
102
103 virtual double GetNoDataValue( int *pbSuccess = NULL );
104
105 void FindPDSTemplate();
106
107 void UncacheData();
108
109 private:
110
111 CPLErr LoadData();
112
113 static void ReadGribData( DataSource &, sInt4, int, double**, grib_MetaData**);
114 sInt4 start;
115 int subgNum;
116 char *longFstLevel;
117
118 double * m_Grib_Data;
119 grib_MetaData* m_Grib_MetaData;
120
121 int nGribDataXSize;
122 int nGribDataYSize;
123 };
124
125 /************************************************************************/
126 /* ConvertUnitInText() */
127 /************************************************************************/
128
ConvertUnitInText(int bMetricUnits,const char * pszTxt)129 static CPLString ConvertUnitInText(int bMetricUnits, const char* pszTxt)
130 {
131 if( !bMetricUnits )
132 return pszTxt;
133
134 CPLString osRes(pszTxt);
135 size_t iPos = osRes.find("[K]");
136 if( iPos != std::string::npos )
137 osRes = osRes.substr(0, iPos) + "[C]" + osRes.substr(iPos + 3);
138 return osRes;
139 }
140
141 /************************************************************************/
142 /* GRIBRasterBand() */
143 /************************************************************************/
144
GRIBRasterBand(GRIBDataset * poDS,int nBand,inventoryType * psInv)145 GRIBRasterBand::GRIBRasterBand( GRIBDataset *poDS, int nBand,
146 inventoryType *psInv )
147 : m_Grib_Data(NULL), m_Grib_MetaData(NULL)
148 {
149 this->poDS = poDS;
150 this->nBand = nBand;
151 this->start = psInv->start;
152 this->subgNum = psInv->subgNum;
153 this->longFstLevel = CPLStrdup(psInv->longFstLevel);
154
155 eDataType = GDT_Float64; // let user do -ot Float32 if needed for saving space, GRIB contains Float64 (though not fully utilized most of the time)
156
157 nBlockXSize = poDS->nRasterXSize;
158 nBlockYSize = 1;
159
160 nGribDataXSize = poDS->nRasterXSize;
161 nGribDataYSize = poDS->nRasterYSize;
162
163 const char* pszGribNormalizeUnits = CPLGetConfigOption("GRIB_NORMALIZE_UNITS", "YES");
164 int bMetricUnits = CSLTestBoolean(pszGribNormalizeUnits);
165
166 SetMetadataItem( "GRIB_UNIT", ConvertUnitInText(bMetricUnits, psInv->unitName) );
167 SetMetadataItem( "GRIB_COMMENT", ConvertUnitInText(bMetricUnits, psInv->comment) );
168 SetMetadataItem( "GRIB_ELEMENT", psInv->element );
169 SetMetadataItem( "GRIB_SHORT_NAME", psInv->shortFstLevel );
170 SetMetadataItem( "GRIB_REF_TIME",
171 CPLString().Printf("%12.0f sec UTC", psInv->refTime ) );
172 SetMetadataItem( "GRIB_VALID_TIME",
173 CPLString().Printf("%12.0f sec UTC", psInv->validTime ) );
174 SetMetadataItem( "GRIB_FORECAST_SECONDS",
175 CPLString().Printf("%.0f sec", psInv->foreSec ) );
176 }
177
178 /************************************************************************/
179 /* FindPDSTemplate() */
180 /* */
181 /* Scan the file for the PDS template info and represent it as */
182 /* metadata. */
183 /************************************************************************/
184
FindPDSTemplate()185 void GRIBRasterBand::FindPDSTemplate()
186
187 {
188 GRIBDataset *poGDS = (GRIBDataset *) poDS;
189
190 /* -------------------------------------------------------------------- */
191 /* Collect section 4 octet information ... we read the file */
192 /* ourselves since the GRIB API does not appear to preserve all */
193 /* this for us. */
194 /* -------------------------------------------------------------------- */
195 GIntBig nOffset = VSIFTellL( poGDS->fp );
196 GByte abyHead[5];
197 GUInt32 nSectSize;
198
199 VSIFSeekL( poGDS->fp, start+16, SEEK_SET );
200 VSIFReadL( abyHead, 5, 1, poGDS->fp );
201
202 while( abyHead[4] != 4 )
203 {
204 memcpy( &nSectSize, abyHead, 4 );
205 CPL_MSBPTR32( &nSectSize );
206
207 if( VSIFSeekL( poGDS->fp, nSectSize-5, SEEK_CUR ) != 0
208 || VSIFReadL( abyHead, 5, 1, poGDS->fp ) != 1 )
209 break;
210 }
211
212 if( abyHead[4] == 4 )
213 {
214 GUInt16 nCoordCount;
215 GUInt16 nPDTN;
216 CPLString osOctet;
217 int i;
218 GByte *pabyBody;
219
220 memcpy( &nSectSize, abyHead, 4 );
221 CPL_MSBPTR32( &nSectSize );
222
223 pabyBody = (GByte *) CPLMalloc(nSectSize-5);
224 VSIFReadL( pabyBody, 1, nSectSize-5, poGDS->fp );
225
226 memcpy( &nCoordCount, pabyBody + 5 - 5, 2 );
227 CPL_MSBPTR16( &nCoordCount );
228
229 memcpy( &nPDTN, pabyBody + 7 - 5, 2 );
230 CPL_MSBPTR16( &nPDTN );
231
232 SetMetadataItem( "GRIB_PDS_PDTN",
233 CPLString().Printf( "%d", nPDTN ) );
234
235 for( i = 9; i < (int) nSectSize; i++ )
236 {
237 char szByte[10];
238
239 if( i == 9 )
240 sprintf( szByte, "%d", pabyBody[i-5] );
241 else
242 sprintf( szByte, " %d", pabyBody[i-5] );
243 osOctet += szByte;
244 }
245
246 SetMetadataItem( "GRIB_PDS_TEMPLATE_NUMBERS", osOctet );
247
248 CPLFree( pabyBody );
249 }
250
251 VSIFSeekL( poGDS->fp, nOffset, SEEK_SET );
252 }
253
254 /************************************************************************/
255 /* GetDescription() */
256 /************************************************************************/
257
GetDescription() const258 const char * GRIBRasterBand::GetDescription() const
259 {
260 if( longFstLevel == NULL )
261 return GDALPamRasterBand::GetDescription();
262 else
263 return longFstLevel;
264 }
265
266 /************************************************************************/
267 /* LoadData() */
268 /************************************************************************/
269
LoadData()270 CPLErr GRIBRasterBand::LoadData()
271
272 {
273 if( !m_Grib_Data )
274 {
275 GRIBDataset *poGDS = (GRIBDataset *) poDS;
276
277 if (poGDS->bCacheOnlyOneBand)
278 {
279 /* In "one-band-at-a-time" strategy, if the last recently used */
280 /* band is not that one, uncache it. We could use a smarter strategy */
281 /* based on a LRU, but that's a bit overkill for now. */
282 poGDS->poLastUsedBand->UncacheData();
283 poGDS->nCachedBytes = 0;
284 }
285 else
286 {
287 /* Once we have cached more than nCachedBytesThreshold bytes, we will switch */
288 /* to "one-band-at-a-time" strategy, instead of caching all bands that have */
289 /* been accessed */
290 if (poGDS->nCachedBytes > poGDS->nCachedBytesThreshold)
291 {
292 CPLDebug("GRIB", "Maximum band cache size reached for this dataset. "
293 "Caching only one band at a time from now");
294 for(int i=0;i<poGDS->nBands;i++)
295 {
296 ((GRIBRasterBand*) poGDS->GetRasterBand(i+1))->UncacheData();
297 }
298 poGDS->nCachedBytes = 0;
299 poGDS->bCacheOnlyOneBand = TRUE;
300 }
301 }
302
303 FileDataSource grib_fp (poGDS->fp);
304
305 // we don't seem to have any way to detect errors in this!
306 ReadGribData(grib_fp, start, subgNum, &m_Grib_Data, &m_Grib_MetaData);
307 if( !m_Grib_Data )
308 {
309 CPLError( CE_Failure, CPLE_AppDefined, "Out of memory." );
310 return CE_Failure;
311 }
312
313 /* -------------------------------------------------------------------- */
314 /* Check that this band matches the dataset as a whole, size */
315 /* wise. (#3246) */
316 /* -------------------------------------------------------------------- */
317 nGribDataXSize = m_Grib_MetaData->gds.Nx;
318 nGribDataYSize = m_Grib_MetaData->gds.Ny;
319
320 poGDS->nCachedBytes += nGribDataXSize * nGribDataYSize * sizeof(double);
321 poGDS->poLastUsedBand = this;
322
323 if( nGribDataXSize != nRasterXSize
324 || nGribDataYSize != nRasterYSize )
325 {
326 CPLError( CE_Warning, CPLE_AppDefined,
327 "Band %d of GRIB dataset is %dx%d, while the first band and dataset is %dx%d. Georeferencing of band %d may be incorrect, and data access may be incomplete.",
328 nBand,
329 nGribDataXSize, nGribDataYSize,
330 nRasterXSize, nRasterYSize,
331 nBand );
332 }
333 }
334
335 return CE_None;
336 }
337
338 /************************************************************************/
339 /* IReadBlock() */
340 /************************************************************************/
341
IReadBlock(CPL_UNUSED int nBlockXOff,int nBlockYOff,void * pImage)342 CPLErr GRIBRasterBand::IReadBlock( CPL_UNUSED int nBlockXOff,
343 int nBlockYOff,
344 void * pImage )
345
346 {
347 CPLErr eErr = LoadData();
348 if (eErr != CE_None)
349 return eErr;
350
351 /* -------------------------------------------------------------------- */
352 /* The image as read is always upside down to our normal */
353 /* orientation so we need to effectively flip it at this */
354 /* point. We also need to deal with bands that are a different */
355 /* size than the dataset as a whole. */
356 /* -------------------------------------------------------------------- */
357 if( nGribDataXSize == nRasterXSize
358 && nGribDataYSize == nRasterYSize )
359 {
360 // Simple 1:1 case.
361 memcpy(pImage,
362 m_Grib_Data + nRasterXSize * (nRasterYSize - nBlockYOff - 1),
363 nRasterXSize * sizeof(double));
364
365 return CE_None;
366 }
367 else
368 {
369 memset( pImage, 0, sizeof(double)*nRasterXSize );
370
371 if( nBlockYOff >= nGribDataYSize ) // off image?
372 return CE_None;
373
374 int nCopyWords = MIN(nRasterXSize,nGribDataXSize);
375
376 memcpy( pImage,
377 m_Grib_Data + nGribDataXSize*(nGribDataYSize-nBlockYOff-1),
378 nCopyWords * sizeof(double) );
379
380 return CE_None;
381 }
382 }
383
384 /************************************************************************/
385 /* GetNoDataValue() */
386 /************************************************************************/
387
GetNoDataValue(int * pbSuccess)388 double GRIBRasterBand::GetNoDataValue( int *pbSuccess )
389 {
390 CPLErr eErr = LoadData();
391 if (eErr != CE_None ||
392 m_Grib_MetaData == NULL ||
393 m_Grib_MetaData->gridAttrib.f_miss == 0)
394 {
395 if (pbSuccess)
396 *pbSuccess = FALSE;
397 return 0;
398 }
399
400 if (m_Grib_MetaData->gridAttrib.f_miss == 2)
401 {
402 /* what TODO ? */
403 CPLDebug("GRIB", "Secondary missing value also set for band %d : %f",
404 nBand, m_Grib_MetaData->gridAttrib.missSec);
405 }
406
407 if (pbSuccess)
408 *pbSuccess = TRUE;
409 return m_Grib_MetaData->gridAttrib.missPri;
410 }
411
412 /************************************************************************/
413 /* ReadGribData() */
414 /************************************************************************/
415
ReadGribData(DataSource & fp,sInt4 start,int subgNum,double ** data,grib_MetaData ** metaData)416 void GRIBRasterBand::ReadGribData( DataSource & fp, sInt4 start, int subgNum, double** data, grib_MetaData** metaData)
417 {
418 /* Initialisation, for calling the ReadGrib2Record function */
419 sInt4 f_endMsg = 1; /* 1 if we read the last grid in a GRIB message, or we haven't read any messages. */
420 // int subgNum = 0; /* The subgrid in the message that we are interested in. */
421 sChar f_unit = 2; /* None = 0, English = 1, Metric = 2 */
422 double majEarth = 0; /* -radEarth if < 6000 ignore, otherwise use this to
423 * override the radEarth in the GRIB1 or GRIB2
424 * message. Needed because NCEP uses 6371.2 but GRIB1 could only state 6367.47. */
425 double minEarth = 0; /* -minEarth if < 6000 ignore, otherwise use this to
426 * override the minEarth in the GRIB1 or GRIB2 message. */
427 sChar f_SimpleVer = 4; /* Which version of the simple NDFD Weather table to
428 * use. (1 is 6/2003) (2 is 1/2004) (3 is 2/2004)
429 * (4 is 11/2004) (default 4) */
430 LatLon lwlf; /* lower left corner (cookie slicing) -lwlf */
431 LatLon uprt; /* upper right corner (cookie slicing) -uprt */
432 IS_dataType is; /* Un-parsed meta data for this GRIB2 message. As well as some memory used by the unpacker. */
433
434 lwlf.lat = -100; // lat == -100 instructs the GRIB decoder that we don't want a subgrid
435
436 IS_Init (&is);
437
438 const char* pszGribNormalizeUnits = CPLGetConfigOption("GRIB_NORMALIZE_UNITS", "YES");
439 if ( !CSLTestBoolean(pszGribNormalizeUnits) )
440 f_unit = 0; /* do not normalize units to metric */
441
442 /* Read GRIB message from file position "start". */
443 fp.DataSourceFseek(start, SEEK_SET);
444 uInt4 grib_DataLen = 0; /* Size of Grib_Data. */
445 *metaData = new grib_MetaData();
446 MetaInit (*metaData);
447 ReadGrib2Record (fp, f_unit, data, &grib_DataLen, *metaData, &is, subgNum,
448 majEarth, minEarth, f_SimpleVer, &f_endMsg, &lwlf, &uprt);
449
450 char * errMsg = errSprintf(NULL); // no intention to show errors, just swallow it and free the memory
451 if( errMsg != NULL )
452 CPLDebug( "GRIB", "%s", errMsg );
453 free(errMsg);
454 IS_Free(&is);
455 }
456
457 /************************************************************************/
458 /* UncacheData() */
459 /************************************************************************/
460
UncacheData()461 void GRIBRasterBand::UncacheData()
462 {
463 if (m_Grib_Data)
464 free (m_Grib_Data);
465 m_Grib_Data = NULL;
466 if (m_Grib_MetaData)
467 {
468 MetaFree( m_Grib_MetaData );
469 delete m_Grib_MetaData;
470 }
471 m_Grib_MetaData = NULL;
472 }
473
474 /************************************************************************/
475 /* ~GRIBRasterBand() */
476 /************************************************************************/
477
~GRIBRasterBand()478 GRIBRasterBand::~GRIBRasterBand()
479 {
480 CPLFree(longFstLevel);
481 UncacheData();
482 }
483
484 /************************************************************************/
485 /* ==================================================================== */
486 /* GRIBDataset */
487 /* ==================================================================== */
488 /************************************************************************/
489
GRIBDataset()490 GRIBDataset::GRIBDataset()
491
492 {
493 poTransform = NULL;
494 pszProjection = CPLStrdup("");
495 adfGeoTransform[0] = 0.0;
496 adfGeoTransform[1] = 1.0;
497 adfGeoTransform[2] = 0.0;
498 adfGeoTransform[3] = 0.0;
499 adfGeoTransform[4] = 0.0;
500 adfGeoTransform[5] = 1.0;
501
502 nCachedBytes = 0;
503 /* Switch caching strategy once 100 MB threshold is reached */
504 /* Why 100 MB ? --> why not ! */
505 nCachedBytesThreshold = ((GIntBig)atoi(CPLGetConfigOption("GRIB_CACHEMAX", "100"))) * 1024 * 1024;
506 bCacheOnlyOneBand = FALSE;
507 poLastUsedBand = NULL;
508 }
509
510 /************************************************************************/
511 /* ~GRIBDataset() */
512 /************************************************************************/
513
~GRIBDataset()514 GRIBDataset::~GRIBDataset()
515
516 {
517 FlushCache();
518 if( fp != NULL )
519 VSIFCloseL( fp );
520
521 CPLFree( pszProjection );
522 }
523
524 /************************************************************************/
525 /* GetGeoTransform() */
526 /************************************************************************/
527
GetGeoTransform(double * padfTransform)528 CPLErr GRIBDataset::GetGeoTransform( double * padfTransform )
529
530 {
531 memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
532 return CE_None;
533 }
534
535 /************************************************************************/
536 /* GetProjectionRef() */
537 /************************************************************************/
538
GetProjectionRef()539 const char *GRIBDataset::GetProjectionRef()
540
541 {
542 return pszProjection;
543 }
544
545 /************************************************************************/
546 /* Identify() */
547 /************************************************************************/
548
Identify(GDALOpenInfo * poOpenInfo)549 int GRIBDataset::Identify( GDALOpenInfo * poOpenInfo )
550 {
551 if (poOpenInfo->nHeaderBytes < 8)
552 return FALSE;
553
554 /* -------------------------------------------------------------------- */
555 /* Does a part of what ReadSECT0() but in a thread-safe way. */
556 /* -------------------------------------------------------------------- */
557 int i;
558 for(i=0;i<poOpenInfo->nHeaderBytes-3;i++)
559 {
560 if (EQUALN((const char*)poOpenInfo->pabyHeader + i, "GRIB", 4) ||
561 EQUALN((const char*)poOpenInfo->pabyHeader + i, "TDLP", 4))
562 return TRUE;
563 }
564
565 return FALSE;
566 }
567
568 /************************************************************************/
569 /* Open() */
570 /************************************************************************/
571
Open(GDALOpenInfo * poOpenInfo)572 GDALDataset *GRIBDataset::Open( GDALOpenInfo * poOpenInfo )
573
574 {
575 if( !Identify(poOpenInfo) )
576 return NULL;
577
578 /* -------------------------------------------------------------------- */
579 /* A fast "probe" on the header that is partially read in memory. */
580 /* -------------------------------------------------------------------- */
581 char *buff = NULL;
582 uInt4 buffLen = 0;
583 sInt4 sect0[SECT0LEN_WORD];
584 uInt4 gribLen;
585 int version;
586 // grib is not thread safe, make sure not to cause problems
587 // for other thread safe formats
588
589 CPLMutexHolderD(&hGRIBMutex);
590 MemoryDataSource mds (poOpenInfo->pabyHeader, poOpenInfo->nHeaderBytes);
591 if (ReadSECT0 (mds, &buff, &buffLen, -1, sect0, &gribLen, &version) < 0) {
592 free (buff);
593 char * errMsg = errSprintf(NULL);
594 if( errMsg != NULL && strstr(errMsg,"Ran out of file") == NULL )
595 CPLDebug( "GRIB", "%s", errMsg );
596 free(errMsg);
597 return NULL;
598 }
599 free(buff);
600
601 /* -------------------------------------------------------------------- */
602 /* Confirm the requested access is supported. */
603 /* -------------------------------------------------------------------- */
604 if( poOpenInfo->eAccess == GA_Update )
605 {
606 CPLError( CE_Failure, CPLE_NotSupported,
607 "The GRIB driver does not support update access to existing"
608 " datasets.\n" );
609 return NULL;
610 }
611 /* -------------------------------------------------------------------- */
612 /* Create a corresponding GDALDataset. */
613 /* -------------------------------------------------------------------- */
614 GRIBDataset *poDS;
615
616 poDS = new GRIBDataset();
617
618 poDS->fp = VSIFOpenL( poOpenInfo->pszFilename, "r" );
619
620 /* Check the return values */
621 if (!poDS->fp) {
622 // we have no FP, so we don't have anywhere to read from
623 char * errMsg = errSprintf(NULL);
624 if( errMsg != NULL )
625 CPLDebug( "GRIB", "%s", errMsg );
626 free(errMsg);
627
628 CPLError( CE_Failure, CPLE_OpenFailed, "Error (%d) opening file %s", errno, poOpenInfo->pszFilename);
629 CPLReleaseMutex(hGRIBMutex); // Release hGRIBMutex otherwise we'll deadlock with GDALDataset own hGRIBMutex
630 delete poDS;
631 CPLAcquireMutex(hGRIBMutex, 1000.0);
632 return NULL;
633 }
634
635 /* -------------------------------------------------------------------- */
636 /* Read the header. */
637 /* -------------------------------------------------------------------- */
638
639 /* -------------------------------------------------------------------- */
640 /* Make an inventory of the GRIB file. */
641 /* The inventory does not contain all the information needed for */
642 /* creating the RasterBands (especially the x and y size), therefore */
643 /* the first GRIB band is also read for some additional metadata. */
644 /* The band-data that is read is stored into the first RasterBand, */
645 /* simply so that the same portion of the file is not read twice. */
646 /* -------------------------------------------------------------------- */
647
648 VSIFSeekL( poDS->fp, 0, SEEK_SET );
649
650 FileDataSource grib_fp (poDS->fp);
651
652 inventoryType *Inv = NULL; /* Contains an GRIB2 message inventory of the file */
653 uInt4 LenInv = 0; /* size of Inv (also # of GRIB2 messages) */
654 int msgNum =0; /* The messageNumber during the inventory. */
655
656 if (GRIB2Inventory (grib_fp, &Inv, &LenInv, 0, &msgNum) <= 0 )
657 {
658 char * errMsg = errSprintf(NULL);
659 if( errMsg != NULL )
660 CPLDebug( "GRIB", "%s", errMsg );
661 free(errMsg);
662
663 CPLError( CE_Failure, CPLE_OpenFailed,
664 "%s is a grib file, but no raster dataset was successfully identified.",
665 poOpenInfo->pszFilename );
666 CPLReleaseMutex(hGRIBMutex); // Release hGRIBMutex otherwise we'll deadlock with GDALDataset own hGRIBMutex
667 delete poDS;
668 CPLAcquireMutex(hGRIBMutex, 1000.0);
669 return NULL;
670 }
671
672 /* -------------------------------------------------------------------- */
673 /* Create band objects. */
674 /* -------------------------------------------------------------------- */
675 GRIBRasterBand *gribBand;
676 for (uInt4 i = 0; i < LenInv; ++i)
677 {
678 uInt4 bandNr = i+1;
679 if (bandNr == 1)
680 {
681 // important: set DataSet extents before creating first RasterBand in it
682 double * data = NULL;
683 grib_MetaData* metaData;
684 GRIBRasterBand::ReadGribData(grib_fp, 0, Inv[i].subgNum, &data, &metaData);
685 if (data == 0 || metaData->gds.Nx < 1 || metaData->gds.Ny < 1)
686 {
687 CPLError( CE_Failure, CPLE_OpenFailed,
688 "%s is a grib file, but no raster dataset was successfully identified.",
689 poOpenInfo->pszFilename );
690 CPLReleaseMutex(hGRIBMutex); // Release hGRIBMutex otherwise we'll deadlock with GDALDataset own hGRIBMutex
691 delete poDS;
692 CPLAcquireMutex(hGRIBMutex, 1000.0);
693 return NULL;
694 }
695
696 poDS->SetGribMetaData(metaData); // set the DataSet's x,y size, georeference and projection from the first GRIB band
697 gribBand = new GRIBRasterBand( poDS, bandNr, Inv+i);
698
699 if( Inv->GribVersion == 2 )
700 gribBand->FindPDSTemplate();
701
702 gribBand->m_Grib_Data = data;
703 gribBand->m_Grib_MetaData = metaData;
704 }
705 else
706 {
707 gribBand = new GRIBRasterBand( poDS, bandNr, Inv+i );
708 if( CSLTestBoolean( CPLGetConfigOption( "GRIB_PDS_ALL_BANDS", "ON" ) ) )
709 {
710 if( Inv->GribVersion == 2 )
711 gribBand->FindPDSTemplate();
712 }
713 }
714 poDS->SetBand( bandNr, gribBand);
715 GRIB2InventoryFree (Inv + i);
716 }
717 free (Inv);
718
719 /* -------------------------------------------------------------------- */
720 /* Initialize any PAM information. */
721 /* -------------------------------------------------------------------- */
722 poDS->SetDescription( poOpenInfo->pszFilename );
723
724 CPLReleaseMutex(hGRIBMutex); // Release hGRIBMutex otherwise we'll deadlock with GDALDataset own hGRIBMutex
725 poDS->TryLoadXML();
726
727 /* -------------------------------------------------------------------- */
728 /* Check for external overviews. */
729 /* -------------------------------------------------------------------- */
730 poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename, poOpenInfo->GetSiblingFiles() );
731 CPLAcquireMutex(hGRIBMutex, 1000.0);
732
733 return( poDS );
734 }
735
736 /************************************************************************/
737 /* SetMetadata() */
738 /************************************************************************/
739
SetGribMetaData(grib_MetaData * meta)740 void GRIBDataset::SetGribMetaData(grib_MetaData* meta)
741 {
742 nRasterXSize = meta->gds.Nx;
743 nRasterYSize = meta->gds.Ny;
744
745 /* -------------------------------------------------------------------- */
746 /* Image projection. */
747 /* -------------------------------------------------------------------- */
748 OGRSpatialReference oSRS;
749
750 switch(meta->gds.projType)
751 {
752 case GS3_LATLON:
753 case GS3_GAUSSIAN_LATLON:
754 // No projection, only latlon system (geographic)
755 break;
756 case GS3_MERCATOR:
757 oSRS.SetMercator(meta->gds.meshLat, meta->gds.orientLon,
758 1.0, 0.0, 0.0);
759 break;
760 case GS3_POLAR:
761 oSRS.SetPS(meta->gds.meshLat, meta->gds.orientLon,
762 meta->gds.scaleLat1,
763 0.0, 0.0);
764 break;
765 case GS3_LAMBERT:
766 oSRS.SetLCC(meta->gds.scaleLat1, meta->gds.scaleLat2,
767 meta->gds.meshLat, meta->gds.orientLon,
768 0.0, 0.0); // set projection
769 break;
770
771
772 case GS3_ORTHOGRAPHIC:
773
774 //oSRS.SetOrthographic(0.0, meta->gds.orientLon,
775 // meta->gds.lon2, meta->gds.lat2);
776 //oSRS.SetGEOS(meta->gds.orientLon, meta->gds.stretchFactor, meta->gds.lon2, meta->gds.lat2);
777 oSRS.SetGEOS( 0, 35785831, 0, 0 ); // hardcoded for now, I don't know yet how to parse the meta->gds section
778 break;
779 case GS3_EQUATOR_EQUIDIST:
780 break;
781 case GS3_AZIMUTH_RANGE:
782 break;
783 }
784
785 /* -------------------------------------------------------------------- */
786 /* Earth model */
787 /* -------------------------------------------------------------------- */
788 double a = meta->gds.majEarth * 1000.0; // in meters
789 double b = meta->gds.minEarth * 1000.0;
790 if( a == 0 && b == 0 )
791 {
792 a = 6377563.396;
793 b = 6356256.910;
794 }
795
796 if (meta->gds.f_sphere)
797 {
798 oSRS.SetGeogCS( "Coordinate System imported from GRIB file",
799 NULL,
800 "Sphere",
801 a, 0.0 );
802 }
803 else
804 {
805 double fInv = a/(a-b);
806 oSRS.SetGeogCS( "Coordinate System imported from GRIB file",
807 NULL,
808 "Spheroid imported from GRIB file",
809 a, fInv );
810 }
811
812 OGRSpatialReference oLL; // construct the "geographic" part of oSRS
813 oLL.CopyGeogCSFrom( &oSRS );
814
815 double rMinX;
816 double rMaxY;
817 double rPixelSizeX;
818 double rPixelSizeY;
819 if (meta->gds.projType == GS3_ORTHOGRAPHIC)
820 {
821 //rMinX = -meta->gds.Dx * (meta->gds.Nx / 2); // This is what should work, but it doesn't .. Dx seems to have an inverse relation with pixel size
822 //rMaxY = meta->gds.Dy * (meta->gds.Ny / 2);
823 const double geosExtentInMeters = 11137496.552; // hardcoded for now, assumption: GEOS projection, full disc (like MSG)
824 rMinX = -(geosExtentInMeters / 2);
825 rMaxY = geosExtentInMeters / 2;
826 rPixelSizeX = geosExtentInMeters / meta->gds.Nx;
827 rPixelSizeY = geosExtentInMeters / meta->gds.Ny;
828 }
829 else if( oSRS.IsProjected() )
830 {
831 rMinX = meta->gds.lon1; // longitude in degrees, to be transformed to meters (or degrees in case of latlon)
832 rMaxY = meta->gds.lat1; // latitude in degrees, to be transformed to meters
833 OGRCoordinateTransformation *poTransformLLtoSRS = OGRCreateCoordinateTransformation( &(oLL), &(oSRS) );
834 if ((poTransformLLtoSRS != NULL) && poTransformLLtoSRS->Transform( 1, &rMinX, &rMaxY )) // transform it to meters
835 {
836 if (meta->gds.scan == GRIB2BIT_2) // Y is minY, GDAL wants maxY
837 rMaxY += (meta->gds.Ny - 1) * meta->gds.Dy; // -1 because we GDAL needs the coordinates of the centre of the pixel
838 rPixelSizeX = meta->gds.Dx;
839 rPixelSizeY = meta->gds.Dy;
840 }
841 else
842 {
843 rMinX = 0.0;
844 rMaxY = 0.0;
845
846 rPixelSizeX = 1.0;
847 rPixelSizeY = -1.0;
848
849 oSRS.Clear();
850
851 CPLError( CE_Warning, CPLE_AppDefined,
852 "Unable to perform coordinate transformations, so the correct\n"
853 "projected geotransform could not be deduced from the lat/long\n"
854 "control points. Defaulting to ungeoreferenced." );
855 }
856 delete poTransformLLtoSRS;
857 }
858 else
859 {
860 rMinX = meta->gds.lon1; // longitude in degrees, to be transformed to meters (or degrees in case of latlon)
861 rMaxY = meta->gds.lat1; // latitude in degrees, to be transformed to meters
862
863 double rMinY = meta->gds.lat2;
864 if (meta->gds.lat2 > rMaxY)
865 {
866 rMaxY = meta->gds.lat2;
867 rMinY = meta->gds.lat1;
868 }
869
870 if( meta->gds.Nx == 1 )
871 rPixelSizeX = meta->gds.Dx;
872 else if (meta->gds.lon1 > meta->gds.lon2)
873 rPixelSizeX = (360.0 - (meta->gds.lon1 - meta->gds.lon2)) / (meta->gds.Nx - 1);
874 else
875 rPixelSizeX = (meta->gds.lon2 - meta->gds.lon1) / (meta->gds.Nx - 1);
876
877 if( meta->gds.Ny == 1 )
878 rPixelSizeY = meta->gds.Dy;
879 else
880 rPixelSizeY = (rMaxY - rMinY) / (meta->gds.Ny - 1);
881
882 // Do some sanity checks for cases that can't be handled by the above
883 // pixel size corrections. GRIB1 has a minimum precision of 0.001
884 // for latitudes and longitudes, so we'll allow a bit higher than that.
885 if (rPixelSizeX < 0 || fabs(rPixelSizeX - meta->gds.Dx) > 0.002)
886 rPixelSizeX = meta->gds.Dx;
887
888 if (rPixelSizeY < 0 || fabs(rPixelSizeY - meta->gds.Dy) > 0.002)
889 rPixelSizeY = meta->gds.Dy;
890 }
891
892 // http://gdal.org/gdal_datamodel.html :
893 // we need the top left corner of the top left pixel.
894 // At the moment we have the center of the pixel.
895 rMinX-=rPixelSizeX/2;
896 rMaxY+=rPixelSizeY/2;
897
898 adfGeoTransform[0] = rMinX;
899 adfGeoTransform[3] = rMaxY;
900 adfGeoTransform[1] = rPixelSizeX;
901 adfGeoTransform[5] = -rPixelSizeY;
902
903 CPLFree( pszProjection );
904 pszProjection = NULL;
905 oSRS.exportToWkt( &(pszProjection) );
906 }
907
908 /************************************************************************/
909 /* GDALDeregister_GRIB() */
910 /************************************************************************/
911
GDALDeregister_GRIB(GDALDriver *)912 static void GDALDeregister_GRIB(GDALDriver* )
913 {
914 if( hGRIBMutex != NULL )
915 {
916 CPLDestroyMutex(hGRIBMutex);
917 hGRIBMutex = NULL;
918 }
919 }
920
921 /************************************************************************/
922 /* GDALRegister_GRIB() */
923 /************************************************************************/
924
GDALRegister_GRIB()925 void GDALRegister_GRIB()
926
927 {
928 GDALDriver *poDriver;
929
930 if( GDALGetDriverByName( "GRIB" ) == NULL )
931 {
932 poDriver = new GDALDriver();
933
934 poDriver->SetDescription( "GRIB" );
935 poDriver->SetMetadataItem( GDAL_DCAP_RASTER, "YES" );
936 poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
937 "GRIdded Binary (.grb)" );
938 poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
939 "frmt_grib.html" );
940 poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "grb" );
941 poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
942
943 poDriver->pfnOpen = GRIBDataset::Open;
944 poDriver->pfnIdentify = GRIBDataset::Identify;
945 poDriver->pfnUnloadDriver = GDALDeregister_GRIB;
946
947 GetGDALDriverManager()->RegisterDriver( poDriver );
948 }
949 }
950