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