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