1 : /******************************************************************************
2 : * $Id: gribdataset.cpp 18099 2009-11-25 20:02:16Z warmerdam $
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 18099 2009-11-25 20:02:16Z warmerdam $");
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 : FILE *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 :
250 : /* -------------------------------------------------------------------- */
251 : /* Check that this band matches the dataset as a whole, size */
252 : /* wise. (#3246) */
253 : /* -------------------------------------------------------------------- */
254 3 : nGribDataXSize = m_Grib_MetaData->gds.Nx;
255 3 : nGribDataYSize = m_Grib_MetaData->gds.Ny;
256 :
257 3 : if( nGribDataXSize != nRasterXSize
258 : || nGribDataYSize != nRasterYSize )
259 : {
260 : CPLError( CE_Warning, CPLE_AppDefined,
261 : "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.",
262 : nBand,
263 : nGribDataXSize, nGribDataYSize,
264 : nRasterXSize, nRasterYSize,
265 1 : nBand );
266 3 : }
267 : }
268 :
269 : /* -------------------------------------------------------------------- */
270 : /* The image as read is always upside down to our normal */
271 : /* orientation so we need to effectively flip it at this */
272 : /* point. We also need to deal with bands that are a different */
273 : /* size than the dataset as a whole. */
274 : /* -------------------------------------------------------------------- */
275 281 : if( nGribDataXSize == nRasterXSize
276 : && nGribDataYSize == nRasterYSize )
277 : {
278 : // Simple 1:1 case.
279 : memcpy(pImage,
280 : m_Grib_Data + nRasterXSize * (nRasterYSize - nBlockYOff - 1),
281 203 : nRasterXSize * sizeof(double));
282 :
283 203 : return CE_None;
284 : }
285 : else
286 : {
287 78 : memset( pImage, 0, sizeof(double)*nRasterXSize );
288 :
289 78 : if( nBlockYOff >= nGribDataYSize ) // off image?
290 57 : return CE_None;
291 :
292 21 : int nCopyWords = MIN(nRasterXSize,nGribDataXSize);
293 :
294 : memcpy( pImage,
295 : m_Grib_Data + nGribDataXSize*(nGribDataYSize-nBlockYOff-1),
296 21 : nCopyWords * sizeof(double) );
297 :
298 21 : return CE_None;
299 : }
300 : }
301 :
302 : /************************************************************************/
303 : /* ReadGribData() */
304 : /************************************************************************/
305 :
306 6 : void GRIBRasterBand::ReadGribData( DataSource & fp, sInt4 start, int subgNum, double** data, grib_MetaData** metaData)
307 : {
308 : /* Initialisation, for calling the ReadGrib2Record function */
309 6 : sInt4 f_endMsg = 1; /* 1 if we read the last grid in a GRIB message, or we haven't read any messages. */
310 : // int subgNum = 0; /* The subgrid in the message that we are interested in. */
311 6 : sChar f_unit = 2; /* None = 0, English = 1, Metric = 2 */
312 6 : double majEarth = 0; /* -radEarth if < 6000 ignore, otherwise use this to
313 : * override the radEarth in the GRIB1 or GRIB2
314 : * message. Needed because NCEP uses 6371.2 but GRIB1 could only state 6367.47. */
315 6 : double minEarth = 0; /* -minEarth if < 6000 ignore, otherwise use this to
316 : * override the minEarth in the GRIB1 or GRIB2 message. */
317 6 : sChar f_SimpleVer = 4; /* Which version of the simple NDFD Weather table to
318 : * use. (1 is 6/2003) (2 is 1/2004) (3 is 2/2004)
319 : * (4 is 11/2004) (default 4) */
320 : LatLon lwlf; /* lower left corner (cookie slicing) -lwlf */
321 : LatLon uprt; /* upper right corner (cookie slicing) -uprt */
322 : IS_dataType is; /* Un-parsed meta data for this GRIB2 message. As well as some memory used by the unpacker. */
323 :
324 6 : lwlf.lat = -100; // lat == -100 instructs the GRIB decoder that we don't want a subgrid
325 :
326 6 : IS_Init (&is);
327 :
328 :
329 : /* Read GRIB message from file position "start". */
330 6 : fp.DataSourceFseek(start, SEEK_SET);
331 6 : uInt4 grib_DataLen = 0; /* Size of Grib_Data. */
332 6 : *metaData = new grib_MetaData();
333 6 : MetaInit (*metaData);
334 : ReadGrib2Record (fp, f_unit, data, &grib_DataLen, *metaData, &is, subgNum,
335 6 : majEarth, minEarth, f_SimpleVer, &f_endMsg, &lwlf, &uprt);
336 :
337 6 : char * errMsg = errSprintf(NULL); // no intention to show errors, just swallow it and free the memory
338 6 : if( errMsg != NULL )
339 0 : CPLDebug( "GRIB", "%s", errMsg );
340 6 : free(errMsg);
341 6 : IS_Free(&is);
342 6 : }
343 :
344 : /************************************************************************/
345 : /* ~GRIBRasterBand() */
346 : /************************************************************************/
347 :
348 36 : GRIBRasterBand::~GRIBRasterBand()
349 : {
350 18 : CPLFree(longFstLevel);
351 18 : if (m_Grib_Data)
352 6 : free (m_Grib_Data);
353 18 : if (m_Grib_MetaData)
354 : {
355 6 : MetaFree( m_Grib_MetaData );
356 6 : delete m_Grib_MetaData;
357 : }
358 36 : }
359 :
360 : /************************************************************************/
361 : /* ==================================================================== */
362 : /* GRIBDataset */
363 : /* ==================================================================== */
364 : /************************************************************************/
365 :
366 3 : GRIBDataset::GRIBDataset()
367 :
368 : {
369 3 : poTransform = NULL;
370 3 : pszProjection = CPLStrdup("");
371 3 : adfGeoTransform[0] = 0.0;
372 3 : adfGeoTransform[1] = 1.0;
373 3 : adfGeoTransform[2] = 0.0;
374 3 : adfGeoTransform[3] = 0.0;
375 3 : adfGeoTransform[4] = 0.0;
376 3 : adfGeoTransform[5] = 1.0;
377 3 : }
378 :
379 : /************************************************************************/
380 : /* ~GRIBDataset() */
381 : /************************************************************************/
382 :
383 6 : GRIBDataset::~GRIBDataset()
384 :
385 : {
386 3 : FlushCache();
387 3 : if( fp != NULL )
388 3 : VSIFCloseL( fp );
389 :
390 3 : CPLFree( pszProjection );
391 6 : }
392 :
393 : /************************************************************************/
394 : /* GetGeoTransform() */
395 : /************************************************************************/
396 :
397 0 : CPLErr GRIBDataset::GetGeoTransform( double * padfTransform )
398 :
399 : {
400 0 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
401 0 : return CE_None;
402 : }
403 :
404 : /************************************************************************/
405 : /* GetProjectionRef() */
406 : /************************************************************************/
407 :
408 0 : const char *GRIBDataset::GetProjectionRef()
409 :
410 : {
411 0 : return pszProjection;
412 : }
413 :
414 : /************************************************************************/
415 : /* Identify() */
416 : /************************************************************************/
417 :
418 8915 : int GRIBDataset::Identify( GDALOpenInfo * poOpenInfo )
419 : {
420 8915 : if (poOpenInfo->nHeaderBytes < 8)
421 8391 : return FALSE;
422 :
423 : /* -------------------------------------------------------------------- */
424 : /* Does a part of what ReadSECT0() but in a thread-safe way. */
425 : /* -------------------------------------------------------------------- */
426 : int i;
427 401626 : for(i=0;i<poOpenInfo->nHeaderBytes-3;i++)
428 : {
429 401105 : if (EQUALN((const char*)poOpenInfo->pabyHeader + i, "GRIB", 4) ||
430 : EQUALN((const char*)poOpenInfo->pabyHeader + i, "TDLP", 4))
431 3 : return TRUE;
432 : }
433 :
434 521 : return FALSE;
435 : }
436 :
437 : /************************************************************************/
438 : /* Open() */
439 : /************************************************************************/
440 :
441 1188 : GDALDataset *GRIBDataset::Open( GDALOpenInfo * poOpenInfo )
442 :
443 : {
444 1188 : if( !Identify(poOpenInfo) )
445 1185 : return NULL;
446 :
447 : /* -------------------------------------------------------------------- */
448 : /* A fast "probe" on the header that is partially read in memory. */
449 : /* -------------------------------------------------------------------- */
450 3 : if( poOpenInfo->fp == NULL)
451 0 : return NULL;
452 :
453 3 : char *buff = NULL;
454 3 : uInt4 buffLen = 0;
455 : sInt4 sect0[SECT0LEN_WORD];
456 : uInt4 gribLen;
457 : int version;
458 : // grib is not thread safe, make sure not to cause problems
459 : // for other thread safe formats
460 : static void *mutex = 0;
461 3 : CPLMutexHolderD(&mutex);
462 3 : MemoryDataSource mds (poOpenInfo->pabyHeader, poOpenInfo->nHeaderBytes);
463 3 : if (ReadSECT0 (mds, &buff, &buffLen, -1, sect0, &gribLen, &version) < 0) {
464 0 : free (buff);
465 0 : char * errMsg = errSprintf(NULL);
466 0 : if( errMsg != NULL && strstr(errMsg,"Ran out of file") == NULL )
467 0 : CPLDebug( "GRIB", "%s", errMsg );
468 0 : free(errMsg);
469 0 : return NULL;
470 : }
471 3 : free(buff);
472 :
473 : /* -------------------------------------------------------------------- */
474 : /* Confirm the requested access is supported. */
475 : /* -------------------------------------------------------------------- */
476 3 : if( poOpenInfo->eAccess == GA_Update )
477 : {
478 : CPLError( CE_Failure, CPLE_NotSupported,
479 : "The GRIB driver does not support update access to existing"
480 0 : " datasets.\n" );
481 0 : return NULL;
482 : }
483 : /* -------------------------------------------------------------------- */
484 : /* Create a corresponding GDALDataset. */
485 : /* -------------------------------------------------------------------- */
486 : GRIBDataset *poDS;
487 :
488 3 : poDS = new GRIBDataset();
489 :
490 6 : poDS->fp = VSIFOpenL( poOpenInfo->pszFilename, "r" );
491 :
492 : /* -------------------------------------------------------------------- */
493 : /* Read the header. */
494 : /* -------------------------------------------------------------------- */
495 :
496 : /* -------------------------------------------------------------------- */
497 : /* Make an inventory of the GRIB file. */
498 : /* The inventory does not contain all the information needed for */
499 : /* creating the RasterBands (especially the x and y size), therefore */
500 : /* the first GRIB band is also read for some additional metadata. */
501 : /* The band-data that is read is stored into the first RasterBand, */
502 : /* simply so that the same portion of the file is not read twice. */
503 : /* -------------------------------------------------------------------- */
504 :
505 3 : VSIFSeekL( poDS->fp, 0, SEEK_SET );
506 :
507 3 : FileDataSource grib_fp (poDS->fp);
508 :
509 3 : inventoryType *Inv = NULL; /* Contains an GRIB2 message inventory of the file */
510 3 : uInt4 LenInv = 0; /* size of Inv (also # of GRIB2 messages) */
511 3 : int msgNum =0; /* The messageNumber during the inventory. */
512 :
513 3 : if (GRIB2Inventory (grib_fp, &Inv, &LenInv, 0, &msgNum) <= 0 )
514 : {
515 0 : char * errMsg = errSprintf(NULL);
516 0 : if( errMsg != NULL )
517 0 : CPLDebug( "GRIB", "%s", errMsg );
518 0 : free(errMsg);
519 :
520 : CPLError( CE_Failure, CPLE_OpenFailed,
521 : "%s is a grib file, but no raster dataset was successfully identified.",
522 0 : poOpenInfo->pszFilename );
523 0 : delete poDS;
524 0 : return NULL;
525 : }
526 :
527 : /* -------------------------------------------------------------------- */
528 : /* Create band objects. */
529 : /* -------------------------------------------------------------------- */
530 21 : for (uInt4 i = 0; i < LenInv; ++i)
531 : {
532 18 : uInt4 bandNr = i+1;
533 18 : if (bandNr == 1)
534 : {
535 : // important: set DataSet extents before creating first RasterBand in it
536 3 : double * data = NULL;
537 : grib_MetaData* metaData;
538 3 : GRIBRasterBand::ReadGribData(grib_fp, 0, Inv[i].subgNum, &data, &metaData);
539 3 : if (metaData->gds.Nx < 1 || metaData->gds.Ny < 1 )
540 : {
541 : CPLError( CE_Failure, CPLE_OpenFailed,
542 : "%s is a grib file, but no raster dataset was successfully identified.",
543 0 : poOpenInfo->pszFilename );
544 0 : delete poDS;
545 0 : return NULL;
546 : }
547 :
548 3 : poDS->SetGribMetaData(metaData); // set the DataSet's x,y size, georeference and projection from the first GRIB band
549 3 : GRIBRasterBand* gribBand = new GRIBRasterBand( poDS, bandNr, Inv+i);
550 :
551 3 : if( Inv->GribVersion == 2 )
552 1 : gribBand->FindPDSTemplate();
553 :
554 3 : gribBand->m_Grib_Data = data;
555 3 : gribBand->m_Grib_MetaData = metaData;
556 3 : poDS->SetBand( bandNr, gribBand);
557 : }
558 : else
559 15 : poDS->SetBand( bandNr, new GRIBRasterBand( poDS, bandNr, Inv+i ));
560 18 : GRIB2InventoryFree (Inv + i);
561 : }
562 3 : free (Inv);
563 :
564 : /* -------------------------------------------------------------------- */
565 : /* Initialize any PAM information. */
566 : /* -------------------------------------------------------------------- */
567 3 : poDS->SetDescription( poOpenInfo->pszFilename );
568 3 : poDS->TryLoadXML();
569 :
570 3 : return( poDS );
571 : }
572 :
573 : /************************************************************************/
574 : /* SetMetadata() */
575 : /************************************************************************/
576 :
577 3 : void GRIBDataset::SetGribMetaData(grib_MetaData* meta)
578 : {
579 3 : nRasterXSize = meta->gds.Nx;
580 3 : nRasterYSize = meta->gds.Ny;
581 :
582 : /* -------------------------------------------------------------------- */
583 : /* Image projection. */
584 : /* -------------------------------------------------------------------- */
585 3 : OGRSpatialReference oSRS;
586 :
587 3 : switch(meta->gds.projType)
588 : {
589 : case GS3_LATLON:
590 : case GS3_GAUSSIAN_LATLON:
591 : // No projection, only latlon system (geographic)
592 2 : break;
593 : case GS3_MERCATOR:
594 : oSRS.SetMercator(meta->gds.meshLat, meta->gds.orientLon,
595 1 : 1.0, 0.0, 0.0);
596 1 : break;
597 : case GS3_POLAR:
598 : oSRS.SetPS(meta->gds.meshLat, meta->gds.orientLon,
599 : meta->gds.scaleLat1,
600 0 : 0.0, 0.0);
601 0 : break;
602 : case GS3_LAMBERT:
603 : oSRS.SetLCC(meta->gds.scaleLat1, meta->gds.scaleLat2,
604 : 0.0, meta->gds.orientLon,
605 0 : 0.0, 0.0); // set projection
606 0 : break;
607 :
608 :
609 : case GS3_ORTHOGRAPHIC:
610 :
611 : //oSRS.SetOrthographic(0.0, meta->gds.orientLon,
612 : // meta->gds.lon2, meta->gds.lat2);
613 : //oSRS.SetGEOS(meta->gds.orientLon, meta->gds.stretchFactor, meta->gds.lon2, meta->gds.lat2);
614 0 : oSRS.SetGEOS( 0, 35785831, 0, 0 ); // hardcoded for now, I don't know yet how to parse the meta->gds section
615 : break;
616 : case GS3_EQUATOR_EQUIDIST:
617 : break;
618 : case GS3_AZIMUTH_RANGE:
619 : break;
620 : }
621 :
622 : /* -------------------------------------------------------------------- */
623 : /* Earth model */
624 : /* -------------------------------------------------------------------- */
625 3 : double a = meta->gds.majEarth * 1000.0; // in meters
626 3 : double b = meta->gds.minEarth * 1000.0;
627 3 : if( a == 0 && b == 0 )
628 : {
629 0 : a = 6377563.396;
630 0 : b = 6356256.910;
631 : }
632 :
633 3 : if (meta->gds.f_sphere)
634 : {
635 : oSRS.SetGeogCS( "Coordinate System imported from GRIB file",
636 : NULL,
637 : "Sphere",
638 3 : a, 0.0 );
639 : }
640 : else
641 : {
642 0 : double fInv = a/(a-b);
643 : oSRS.SetGeogCS( "Coordinate System imported from GRIB file",
644 : NULL,
645 : "Spheroid imported from GRIB file",
646 0 : a, fInv );
647 : }
648 :
649 3 : OGRSpatialReference oLL; // construct the "geographic" part of oSRS
650 3 : oLL.CopyGeogCSFrom( &oSRS );
651 :
652 : double rMinX;
653 : double rMaxY;
654 : double rPixelSizeX;
655 : double rPixelSizeY;
656 3 : if (meta->gds.projType == GS3_ORTHOGRAPHIC)
657 : {
658 : //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
659 : //rMaxY = meta->gds.Dy * (meta->gds.Ny / 2);
660 0 : const double geosExtentInMeters = 11137496.552; // hardcoded for now, assumption: GEOS projection, full disc (like MSG)
661 0 : rMinX = -(geosExtentInMeters / 2);
662 0 : rMaxY = geosExtentInMeters / 2;
663 0 : rPixelSizeX = geosExtentInMeters / meta->gds.Nx;
664 0 : rPixelSizeY = geosExtentInMeters / meta->gds.Ny;
665 : }
666 3 : else if( oSRS.IsProjected() )
667 : {
668 1 : rMinX = meta->gds.lon1; // longitude in degrees, to be transformed to meters (or degrees in case of latlon)
669 1 : rMaxY = meta->gds.lat1; // latitude in degrees, to be transformed to meters
670 1 : OGRCoordinateTransformation *poTransformLLtoSRS = OGRCreateCoordinateTransformation( &(oLL), &(oSRS) );
671 1 : if ((poTransformLLtoSRS != NULL) && poTransformLLtoSRS->Transform( 1, &rMinX, &rMaxY )) // transform it to meters
672 : {
673 1 : if (meta->gds.scan == GRIB2BIT_2) // Y is minY, GDAL wants maxY
674 1 : rMaxY += (meta->gds.Ny - 1) * meta->gds.Dy; // -1 because we GDAL needs the coordinates of the centre of the pixel
675 1 : rPixelSizeX = meta->gds.Dx;
676 1 : rPixelSizeY = meta->gds.Dy;
677 : }
678 : else
679 : {
680 0 : rMinX = 0.0;
681 0 : rMaxY = 0.0;
682 :
683 0 : rPixelSizeX = 1.0;
684 0 : rPixelSizeY = -1.0;
685 :
686 0 : oSRS.Clear();
687 :
688 : CPLError( CE_Warning, CPLE_AppDefined,
689 : "Unable to perform coordinate transformations, so the correct\n"
690 : "projected geotransform could not be deduced from the lat/long\n"
691 0 : "control points. Defaulting to ungeoreferenced." );
692 : }
693 1 : delete poTransformLLtoSRS;
694 : }
695 : else
696 : {
697 2 : rMinX = meta->gds.lon1; // longitude in degrees, to be transformed to meters (or degrees in case of latlon)
698 2 : rMaxY = meta->gds.lat1; // latitude in degrees, to be transformed to meters
699 :
700 2 : if (meta->gds.scan == GRIB2BIT_2) // Y is minY, GDAL wants maxY
701 2 : rMaxY += (meta->gds.Ny - 1) * meta->gds.Dy; // -1 because we GDAL needs the coordinates of the centre of the pixel
702 2 : rPixelSizeX = meta->gds.Dx;
703 2 : rPixelSizeY = meta->gds.Dy;
704 : }
705 :
706 3 : adfGeoTransform[0] = rMinX;
707 3 : adfGeoTransform[3] = rMaxY;
708 3 : adfGeoTransform[1] = rPixelSizeX;
709 3 : adfGeoTransform[5] = -rPixelSizeY;
710 :
711 3 : CPLFree( pszProjection );
712 3 : pszProjection = NULL;
713 3 : oSRS.exportToWkt( &(pszProjection) );
714 3 : }
715 :
716 : /************************************************************************/
717 : /* GDALRegister_GRIB() */
718 : /************************************************************************/
719 :
720 338 : void GDALRegister_GRIB()
721 :
722 : {
723 : GDALDriver *poDriver;
724 :
725 338 : if( GDALGetDriverByName( "GRIB" ) == NULL )
726 : {
727 336 : poDriver = new GDALDriver();
728 :
729 336 : poDriver->SetDescription( "GRIB" );
730 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
731 336 : "GRIdded Binary (.grb)" );
732 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
733 336 : "frmt_grib.html" );
734 336 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "grb" );
735 :
736 336 : poDriver->pfnOpen = GRIBDataset::Open;
737 336 : poDriver->pfnIdentify = GRIBDataset::Identify;
738 :
739 336 : GetGDALDriverManager()->RegisterDriver( poDriver );
740 : }
741 338 : }
|
