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