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