Changeset 13823

Show
Ignore:
Timestamp:
02/19/08 16:56:06 (5 months ago)
Author:
warmerdam
Message:

expand tabs, reformat to gdal-ish style

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • spike/grib/gribdataset.cpp

    r13662 r13823  
    8585    friend class GRIBDataset; 
    8686     
    87   public: 
    88         GRIBRasterBand( GRIBDataset*, int, sInt4, int, char* ); 
    89         virtual ~GRIBRasterBand(); 
    90         virtual CPLErr IReadBlock( int, int, void * ); 
    91         virtual const char *GetDescription() const; 
    92        private: 
    93         static void ReadGribData( DataSource &, sInt4, int, double**, grib_MetaData**); 
    94         sInt4 start; 
    95                                int subgNum; 
    96                                char *longFstLevel; 
    97         double * m_Grib_Data; 
    98         grib_MetaData* m_Grib_MetaData; 
     87public: 
     88    GRIBRasterBand( GRIBDataset*, int, sInt4, int, char* ); 
     89    virtual ~GRIBRasterBand(); 
     90    virtual CPLErr IReadBlock( int, int, void * ); 
     91    virtual const char *GetDescription() const; 
     92private: 
     93    static void ReadGribData( DataSource &, sInt4, int, double**, grib_MetaData**); 
     94    sInt4 start; 
     95    int subgNum; 
     96    char *longFstLevel; 
     97    double * m_Grib_Data; 
     98    grib_MetaData* m_Grib_MetaData; 
    9999}; 
    100100 
    …… 
    105105 
    106106GRIBRasterBand::GRIBRasterBand( GRIBDataset *poDS, int nBand, sInt4 start, int subgNum, char *longFstLevel ) 
    107 : m_Grib_Data(NULL) 
    108 , m_Grib_MetaData(NULL) 
     107  : m_Grib_Data(NULL) 
     108  , m_Grib_MetaData(NULL) 
    109109{ 
    110110    this->poDS = poDS; 
    111111    this->nBand = nBand; 
    112                this->start = start; 
    113                this->subgNum = subgNum; 
    114                this->longFstLevel = CPLStrdup(longFstLevel); 
    115  
    116                eDataType = GDT_Float64; // let user do -ot Float32 if needed for saving space, GRIB contains Float64 (though not fully utilized most of the time) 
     112    this->start = start; 
     113    this->subgNum = subgNum; 
     114    this->longFstLevel = CPLStrdup(longFstLevel); 
     115 
     116    eDataType = GDT_Float64; // let user do -ot Float32 if needed for saving space, GRIB contains Float64 (though not fully utilized most of the time) 
    117117 
    118118    nBlockXSize = poDS->nRasterXSize; 
    119119    nBlockYSize = 1; 
    120120 
    121                // printf("Band %d: %s\n", nBand, GetDescription()); 
    122                //printf("size x=%d,y=%d pixels\n", poDS->nRasterXSize, poDS->nRasterYSize); 
    123                //if (m_meta->gridAttrib.f_maxmin) 
    124                //     printf("Min = %f\nMax = %f\n", m_meta->gridAttrib.min, m_meta->gridAttrib.max); 
    125                //if (m_meta->gridAttrib.f_miss) 
    126                //{ 
    127                //     printf("Missing data value = %f\n", m_meta->gridAttrib.missPri); 
    128                //     if (m_meta->gridAttrib.f_miss > 1) 
    129                //             printf("Secondary missing data value = %f\n", m_meta->gridAttrib.missSec); 
    130                //} 
    131                //printf("Content: %s\n", m_meta->longFstLevel); 
    132  
    133                //if (m_meta->gds.numPts != poDS->nRasterYSize * poDS->nRasterXSize) 
    134                //     printf("ERROR: numPts != Nx * Ny? (%ld != %ld * %ld)\n", m_meta->gds.numPts, poDS->nRasterXSize, poDS->nRasterYSize); 
    135                //else if ((m_meta->gds.Dx <= 0) || (m_meta->gds.Dy <= 0)) 
    136                //     printf("Projection code requires Dx (%f) > 0 and Dy (%f) > 0\n", m_meta->gds.Dx, m_meta->gds.Dy); 
     121    // printf("Band %d: %s\n", nBand, GetDescription()); 
     122    //printf("size x=%d,y=%d pixels\n", poDS->nRasterXSize, poDS->nRasterYSize); 
     123    //if (m_meta->gridAttrib.f_maxmin) 
     124    //        printf("Min = %f\nMax = %f\n", m_meta->gridAttrib.min, m_meta->gridAttrib.max); 
     125    //if (m_meta->gridAttrib.f_miss) 
     126    //{ 
     127    //        printf("Missing data value = %f\n", m_meta->gridAttrib.missPri); 
     128    //        if (m_meta->gridAttrib.f_miss > 1) 
     129    //                printf("Secondary missing data value = %f\n", m_meta->gridAttrib.missSec); 
     130    //} 
     131    //printf("Content: %s\n", m_meta->longFstLevel); 
     132 
     133    //if (m_meta->gds.numPts != poDS->nRasterYSize * poDS->nRasterXSize) 
     134    //        printf("ERROR: numPts != Nx * Ny? (%ld != %ld * %ld)\n", m_meta->gds.numPts, poDS->nRasterXSize, poDS->nRasterYSize); 
     135    //else if ((m_meta->gds.Dx <= 0) || (m_meta->gds.Dy <= 0)) 
     136    //        printf("Projection code requires Dx (%f) > 0 and Dy (%f) > 0\n", m_meta->gds.Dx, m_meta->gds.Dy); 
    137137} 
    138138 
    …… 
    154154 
    155155CPLErr GRIBRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff, 
    156                                   void * pImage ) 
    157  
    158 { 
    159                if (!m_Grib_Data) 
    160                { 
    161            GRIBDataset *poGDS = (GRIBDataset *) poDS; 
    162  
    163                FileDataSource grib_fp (poGDS->fp); 
    164  
    165                        ReadGribData(grib_fp, start, subgNum, &m_Grib_Data, &m_Grib_MetaData); 
    166                } 
    167  
    168                // Somehow this decoder guarantees us that the image is upside-down (GRIB scan mode 0100). Therefore reverse it again. 
    169                memcpy(pImage, m_Grib_Data + nRasterXSize * (nRasterYSize - nBlockYOff - 1), nRasterXSize * sizeof(double)); 
     156                                  void * pImage ) 
     157 
     158{ 
     159    if (!m_Grib_Data) 
     160    { 
     161        GRIBDataset *poGDS = (GRIBDataset *) poDS; 
     162 
     163        FileDataSource grib_fp (poGDS->fp); 
     164 
     165        ReadGribData(grib_fp, start, subgNum, &m_Grib_Data, &m_Grib_MetaData); 
     166    } 
     167 
     168    // Somehow this decoder guarantees us that the image is upside-down (GRIB scan mode 0100). Therefore reverse it again. 
     169    memcpy(pImage, m_Grib_Data + nRasterXSize * (nRasterYSize - nBlockYOff - 1), nRasterXSize * sizeof(double)); 
    170170 
    171171    return CE_None; 
    …… 
    178178void GRIBRasterBand::ReadGribData( DataSource & fp, sInt4 start, int subgNum, double** data, grib_MetaData** metaData) 
    179179{ 
    180                /* Initialisation, for calling the ReadGrib2Record function */ 
    181                sInt4 f_endMsg = 1;  /* 1 if we read the last grid in a GRIB message, or we haven't read any messages. */ 
    182                // int subgNum = 0;     /* The subgrid in the message that we are interested in. */ 
    183                sChar f_unit = 2;        /* None = 0, English = 1, Metric = 2 */ 
    184                double majEarth = 0;     /* -radEarth if < 6000 ignore, otherwise use this to 
    185                                                                                                 * override the radEarth in the GRIB1 or GRIB2 
    186                                                                                                 * message.  Needed because NCEP uses 6371.2 but GRIB1 could only state 6367.47. */ 
    187                double minEarth = 0;     /* -minEarth if < 6000 ignore, otherwise use this to 
    188                                                                                                 * override the minEarth in the GRIB1 or GRIB2 message. */ 
    189                sChar f_SimpleVer = 4;   /* Which version of the simple NDFD Weather table to 
    190                                                                                                 * use. (1 is 6/2003) (2 is 1/2004) (3 is 2/2004) 
    191                                                                                                 * (4 is 11/2004) (default 4) */ 
    192                LatLon lwlf;         /* lower left corner (cookie slicing) -lwlf */ 
    193                LatLon uprt;         /* upper right corner (cookie slicing) -uprt */ 
    194                IS_dataType is;      /* Un-parsed meta data for this GRIB2 message. As well as some memory used by the unpacker. */ 
    195  
    196                lwlf.lat = -100; // lat == -100 instructs the GRIB decoder that we don't want a subgrid 
    197  
    198                IS_Init (&is); 
    199  
    200                int bandNr = 0; 
    201  
    202  
    203                /* Read GRIB message from file position "start". */ 
    204        fp.DataSourceFseek(start, SEEK_SET); 
    205                uInt4 grib_DataLen = 0;  /* Size of Grib_Data. */ 
    206                *metaData = new grib_MetaData(); 
    207                MetaInit (*metaData); 
    208                ReadGrib2Record (fp, f_unit, data, &grib_DataLen, *metaData, &is, subgNum, 
    209                                                                                                 majEarth, minEarth, f_SimpleVer, &f_endMsg, &lwlf, &uprt); 
    210  
    211                char * errMsg = errSprintf(NULL); // no intention to show errors, just swallow it and free the memory 
    212                free(errMsg); 
    213                IS_Free(&is); 
     180    /* Initialisation, for calling the ReadGrib2Record function */ 
     181    sInt4 f_endMsg = 1;  /* 1 if we read the last grid in a GRIB message, or we haven't read any messages. */ 
     182    // int subgNum = 0;     /* The subgrid in the message that we are interested in. */ 
     183    sChar f_unit = 2;        /* None = 0, English = 1, Metric = 2 */ 
     184    double majEarth = 0;     /* -radEarth if < 6000 ignore, otherwise use this to 
     185                              * override the radEarth in the GRIB1 or GRIB2 
     186                              * message.  Needed because NCEP uses 6371.2 but GRIB1 could only state 6367.47. */ 
     187    double minEarth = 0;     /* -minEarth if < 6000 ignore, otherwise use this to 
     188                              * override the minEarth in the GRIB1 or GRIB2 message. */ 
     189    sChar f_SimpleVer = 4;   /* Which version of the simple NDFD Weather table to 
     190                              * use. (1 is 6/2003) (2 is 1/2004) (3 is 2/2004) 
     191                              * (4 is 11/2004) (default 4) */ 
     192    LatLon lwlf;         /* lower left corner (cookie slicing) -lwlf */ 
     193    LatLon uprt;         /* upper right corner (cookie slicing) -uprt */ 
     194    IS_dataType is;      /* Un-parsed meta data for this GRIB2 message. As well as some memory used by the unpacker. */ 
     195 
     196    lwlf.lat = -100; // lat == -100 instructs the GRIB decoder that we don't want a subgrid 
     197 
     198    IS_Init (&is); 
     199 
     200    int bandNr = 0; 
     201 
     202 
     203    /* Read GRIB message from file position "start". */ 
     204    fp.DataSourceFseek(start, SEEK_SET); 
     205    uInt4 grib_DataLen = 0;  /* Size of Grib_Data. */ 
     206    *metaData = new grib_MetaData(); 
     207    MetaInit (*metaData); 
     208    ReadGrib2Record (fp, f_unit, data, &grib_DataLen, *metaData, &is, subgNum, 
     209                    majEarth, minEarth, f_SimpleVer, &f_endMsg, &lwlf, &uprt); 
     210 
     211    char * errMsg = errSprintf(NULL); // no intention to show errors, just swallow it and free the memory 
     212    free(errMsg); 
     213    IS_Free(&is); 
    214214} 
    215215 
    …… 
    221221{ 
    222222    CPLFree(longFstLevel); 
    223                if (m_Grib_Data) 
    224                        free (m_Grib_Data); 
    225                if (m_Grib_MetaData) 
    226                        free(m_Grib_MetaData); 
     223    if (m_Grib_Data) 
     224        free (m_Grib_Data); 
     225    if (m_Grib_MetaData) 
     226        free(m_Grib_MetaData); 
    227227} 
    228228 
    …… 
    257257        VSIFClose( fp ); 
    258258                 
    259                CPLFree( pszProjection ); 
     259    CPLFree( pszProjection ); 
    260260} 
    261261 
    …… 
    294294        return NULL; 
    295295 
    296                char *buff = NULL; 
    297                uInt4 buffLen = 0; 
    298                sInt4 sect0[SECT0LEN_WORD]; 
    299                uInt4 gribLen; 
    300                int version; 
    301                MemoryDataSource mds (poOpenInfo->pabyHeader, poOpenInfo->nHeaderBytes); 
     296    char *buff = NULL; 
     297    uInt4 buffLen = 0; 
     298    sInt4 sect0[SECT0LEN_WORD]; 
     299    uInt4 gribLen; 
     300    int version; 
     301    MemoryDataSource mds (poOpenInfo->pabyHeader, poOpenInfo->nHeaderBytes); 
    302302    if (ReadSECT0 (mds, &buff, &buffLen, -1, sect0, &gribLen, &version) < 0) { 
    303                                free (buff); 
    304                                char * errMsg = errSprintf(NULL); 
    305                                free(errMsg); 
    306                                return NULL; 
    307     } 
    308                free(buff); 
     303        free (buff); 
     304        char * errMsg = errSprintf(NULL); 
     305        free(errMsg); 
     306        return NULL; 
     307    } 
     308    free(buff); 
    309309     
    310310/* -------------------------------------------------------------------- */ 
    …… 
    331331/* -------------------------------------------------------------------- */ 
    332332     
    333                VSIFSeek( poDS->fp, 0, SEEK_SET ); 
    334  
    335                FileDataSource grib_fp (poDS->fp); 
     333    VSIFSeek( poDS->fp, 0, SEEK_SET ); 
     334 
     335    FileDataSource grib_fp (poDS->fp); 
    336336 
    337337    inventoryType *Inv = NULL;  /* Contains an GRIB2 message inventory of the file */ 
    …… 
    339339    int msgNum =0;          /* The messageNumber during the inventory. */ 
    340340 
    341                if (GRIB2Inventory (grib_fp, &Inv, &LenInv, 0, &msgNum) >= 0) 
    342                { 
     341    if (GRIB2Inventory (grib_fp, &Inv, &LenInv, 0, &msgNum) >= 0) 
     342    { 
    343343        for (uInt4 i = 0; i < LenInv; ++i) 
    344                                { 
    345                                           uInt4 bandNr = i+1; 
    346             if (bandNr == 1) 
    347                                                 { 
    348                                                          // important: set DataSet extents before creating first RasterBand in it 
    349                                                                double * data = NULL; 
    350                                                                grib_MetaData* metaData; 
    351                                                          GRIBRasterBand::ReadGribData(grib_fp, 0, Inv[i].subgNum, &data, &metaData); 
     344        { 
     345            uInt4 bandNr = i+1; 
     346            if (bandNr == 1) 
     347            { 
     348                // important: set DataSet extents before creating first RasterBand in it 
     349                double * data = NULL; 
     350                grib_MetaData* metaData; 
     351                GRIBRasterBand::ReadGribData(grib_fp, 0, Inv[i].subgNum, &data, &metaData); 
    352352                poDS->SetMetaData(metaData); // set the DataSet's x,y size, georeference and projection from the first GRIB band 
    353                                                                GRIBRasterBand* gribBand = new GRIBRasterBand( poDS, bandNr, Inv[i].start, Inv[i].subgNum, Inv[i].longFstLevel); 
    354                                                                gribBand->m_Grib_Data = data; 
    355                                                                gribBand->m_Grib_MetaData = metaData; 
    356                                                                poDS->SetBand( bandNr, gribBand); 
    357                                                 } 
    358                                                 else 
    359                                                                poDS->SetBand( bandNr, new GRIBRasterBand( poDS, bandNr, Inv[i].start, Inv[i].subgNum, Inv[i].longFstLevel)); 
    360             GRIB2InventoryFree (Inv + i); 
    361                                } 
     353                GRIBRasterBand* gribBand = new GRIBRasterBand( poDS, bandNr, Inv[i].start, Inv[i].subgNum, Inv[i].longFstLevel); 
     354                gribBand->m_Grib_Data = data; 
     355                gribBand->m_Grib_MetaData = metaData; 
     356                poDS->SetBand( bandNr, gribBand); 
     357            } 
     358            else 
     359                poDS->SetBand( bandNr, new GRIBRasterBand( poDS, bandNr, Inv[i].start, Inv[i].subgNum, Inv[i].longFstLevel)); 
     360            GRIB2InventoryFree (Inv + i); 
     361        } 
    362362        free (Inv); 
    363                } 
     363    } 
    364364 
    365365/* -------------------------------------------------------------------- */ 
    …… 
    373373 
    374374/************************************************************************/ 
    375 /*                          GDALRegister_GRIB()                          */ 
     375/*                            SetMetadata()                             */ 
    376376/************************************************************************/ 
    377377 
    378378void GRIBDataset::SetMetaData(grib_MetaData* meta) 
    379379{ 
    380                nRasterXSize = meta->gds.Nx; 
    381                nRasterYSize = meta->gds.Ny; 
    382  
    383                OGRSpatialReference oSRS; // the projection of the image 
    384  
    385                switch(meta->gds.projType) 
    386                { 
    387                case GS3_LATLON: 
    388                case GS3_GAUSSIAN_LATLON: 
    389                        // No projection, only latlon system (geographic) 
    390                        break; 
    391                case GS3_MERCATOR: 
    392                        oSRS.SetMercator(meta->gds.meshLat, meta->gds.orientLon, 
    393                        meta->gds.scaleLat1,  
    394                        0.0, 0.0); 
    395                        break; 
    396                case GS3_POLAR: 
    397                        oSRS.SetPS(meta->gds.meshLat, meta->gds.orientLon, 
    398                  meta->gds.scaleLat1, 
    399                  0.0, 0.0); 
    400                        break; 
    401                case GS3_LAMBERT: 
    402                        oSRS.SetLCC(meta->gds.scaleLat1, meta->gds.scaleLat2, 
    403                                                                0.0, meta->gds.orientLon, 
    404                                                                0.0, 0.0); // set projection 
    405                        break; 
     380    nRasterXSize = meta->gds.Nx; 
     381    nRasterYSize = meta->gds.Ny; 
     382 
     383    OGRSpatialReference oSRS; // the projection of the image 
     384 
     385    switch(meta->gds.projType) 
     386    { 
     387      case GS3_LATLON: 
     388      case GS3_GAUSSIAN_LATLON: 
     389          // No projection, only latlon system (geographic) 
     390          break; 
     391      case GS3_MERCATOR: 
     392        oSRS.SetMercator(meta->gds.meshLat, meta->gds.orientLon, 
     393                         meta->gds.scaleLat1,  
     394                         0.0, 0.0); 
     395        break; 
     396      case GS3_POLAR: 
     397        oSRS.SetPS(meta->gds.meshLat, meta->gds.orientLon, 
     398                   meta->gds.scaleLat1, 
     399                   0.0, 0.0); 
     400        break; 
     401      case GS3_LAMBERT: 
     402        oSRS.SetLCC(meta->gds.scaleLat1, meta->gds.scaleLat2, 
     403                    0.0, meta->gds.orientLon, 
     404                    0.0, 0.0); // set projection 
     405        break; 
    406406                         
    407407 
    408                case GS3_ORTHOGRAPHIC: 
    409  
    410                        //oSRS.SetOrthographic(0.0, meta->gds.orientLon, 
    411                        //                                                                                     meta->gds.lon2, meta->gds.lat2); 
    412                        //oSRS.SetGEOS(meta->gds.orientLon, meta->gds.stretchFactor, meta->gds.lon2, meta->gds.lat2); 
    413                        oSRS.SetGEOS(  0, 35785831, 0, 0 ); // hardcoded for now, I don't know yet how to parse the meta->gds section 
    414                        break; 
    415                case GS3_EQUATOR_EQUIDIST: 
    416                        break; 
    417                case GS3_AZIMUTH_RANGE: 
    418                        break; 
    419                } 
    420  
    421                double a = meta->gds.majEarth * 1000.0; // in meters 
    422                double b = meta->gds.minEarth * 1000.0; 
    423                 if( a == 0 && b == 0 ) 
    424                 { 
    425                     a = 6377563.396; 
    426                     b = 6356256.910; 
    427                 } 
    428  
    429                if (meta->gds.f_sphere) 
    430                { 
    431                        oSRS.SetGeogCS( "Coordinate System imported from GRIB file", 
    432                                                                                                                NULL, 
    433                                                                                                                "Sphere", 
    434                                                                                                                a, 0.0, 
    435                                                                                                                "Greenwich", 0.0, 
    436                                                                                                                NULL, 0.0); 
    437                } 
    438                else 
    439                { 
    440                        double fInv = a/(a-b); 
    441                        oSRS.SetGeogCS( "Coordinate System imported from GRIB file", 
    442                                                                                                                NULL, 
    443                                                                                                                "Spheroid imported from GRIB file", 
    444                                                                                                                a, fInv, 
    445                                                                                                                "Greenwich", 0.0, 
    446                                                                                                                NULL, 0.0); 
    447                } 
    448  
    449                OGRSpatialReference oLL; // construct the "geographic" part of oSRS 
    450                char *pszLL; 
    451                (oSRS.GetAttrNode("GEOGCS"))->exportToWkt(&pszLL); 
    452                oLL.importFromWkt(&pszLL); // set ellipsoid only (latlon system) 
    453                // delete pszLL; 
    454  
    455                double rMinX; 
    456                double rMaxY; 
    457                double rPixelSizeX; 
    458                double rPixelSizeY; 
    459                if (meta->gds.projType == GS3_ORTHOGRAPHIC) 
    460                { 
    461                        //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 
    462                        //rMaxY = meta->gds.Dy * (meta->gds.Ny / 2); 
    463                        const double geosExtentInMeters = 11137496.552; // hardcoded for now, assumption: GEOS projection, full disc (like MSG) 
    464                        rMinX = -(geosExtentInMeters / 2); 
    465                        rMaxY = geosExtentInMeters / 2; 
    466                        rPixelSizeX = geosExtentInMeters / meta->gds.Nx; 
    467                        rPixelSizeY = geosExtentInMeters / meta->gds.Ny; 
    468                } 
    469                else 
    470                { 
    471                        rMinX = meta->gds.lon1; // longitude in degrees, to be transformed to meters (or degrees in case of latlon) 
    472                        rMaxY = meta->gds.lat1; // latitude in degrees, to be transformed to meters  
    473                        OGRCoordinateTransformation *poTransformLLtoSRS = OGRCreateCoordinateTransformation( &(oLL), &(oSRS) ); 
    474                        if ((poTransformLLtoSRS != NULL) && poTransformLLtoSRS->Transform( 1, &rMinX, &rMaxY )) // transform it to meters 
    475                        { 
    476                                if (meta->gds.scan == GRIB2BIT_2) // Y is minY, GDAL wants maxY 
    477                                        rMaxY += (meta->gds.Ny - 1) * meta->gds.Dy; // -1 because we GDAL needs the coordinates of the centre of the pixel 
    478                        } 
    479                        delete poTransformLLtoSRS; 
    480                        rPixelSizeX = meta->gds.Dx; 
    481                        rPixelSizeY = meta->gds.Dy; 
    482                } 
    483                adfGeoTransform[0] = rMinX; 
    484                adfGeoTransform[3] = rMaxY; 
    485                adfGeoTransform[1] = rPixelSizeX; 
    486                adfGeoTransform[5] = -rPixelSizeY; 
    487  
    488                oSRS.exportToWkt( &(pszProjection) ); 
    489 } 
    490  
    491 /************************************************************************/ 
    492 /*                          GDALRegister_GRIB()                          */ 
     408      case GS3_ORTHOGRAPHIC: 
     409 
     410        //oSRS.SetOrthographic(0.0, meta->gds.orientLon, 
     411        //                                                                                    meta->gds.lon2, meta->gds.lat2); 
     412        //oSRS.SetGEOS(meta->gds.orientLon, meta->gds.stretchFactor, meta->gds.lon2, meta->gds.lat2); 
     413        oSRS.SetGEOS(  0, 35785831, 0, 0 ); // hardcoded for now, I don't know yet how to parse the meta->gds section 
     414        break; 
     415      case GS3_EQUATOR_EQUIDIST: 
     416        break; 
     417      case GS3_AZIMUTH_RANGE: 
     418        break; 
     419    } 
     420 
     421    double a = meta->gds.majEarth * 1000.0; // in meters 
     422    double b = meta->gds.minEarth * 1000.0; 
     423    if( a == 0 && b == 0 ) 
     424    { 
     425        a = 6377563.396; 
     426        b = 6356256.910; 
     427    } 
     428 
     429    if (meta->gds.f_sphere) 
     430    { 
     431        oSRS.SetGeogCS( "Coordinate System imported from GRIB file", 
     432                        NULL, 
     433                        "Sphere", 
     434                        a, 0.0, 
     435                        "Greenwich", 0.0, 
     436                        NULL, 0.0); 
     437    } 
     438    else 
     439    { 
     440        double fInv = a/(a-b); 
     441        oSRS.SetGeogCS( "Coordinate System imported from GRIB file", 
     442                        NULL, 
     443                        "Spheroid imported from GRIB file", 
     444                        a, fInv, 
     445                        "Greenwich", 0.0, 
     446                        NULL, 0.0); 
     447    } 
     448 
     449    OGRSpatialReference oLL; // construct the "geographic" part of oSRS 
     450    char *pszLL; 
     451    (oSRS.GetAttrNode("GEOGCS"))->exportToWkt(&pszLL); 
     452    oLL.importFromWkt(&pszLL); // set ellipsoid only (latlon system) 
     453    // delete pszLL; 
     454 
     455    double rMinX; 
     456    double rMaxY; 
     457    double rPixelSizeX; 
     458    double rPixelSizeY; 
     459    if (meta->gds.projType == GS3_ORTHOGRAPHIC) 
     460    { 
     461        //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 
     462        //rMaxY = meta->gds.Dy * (meta->gds.Ny / 2); 
     463        const double geosExtentInMeters = 11137496.552; // hardcoded for now, assumption: GEOS projection, full disc (like MSG) 
     464        rMinX = -(geosExtentInMeters / 2); 
     465        rMaxY = geosExtentInMeters / 2; 
     466        rPixelSizeX = geosExtentInMeters / meta->gds.Nx; 
     467        rPixelSizeY = geosExtentInMeters / meta->gds.Ny; 
     468    } 
     469    else 
     470    { 
     471        rMinX = meta->gds.lon1; // longitude in degrees, to be transformed to meters (or degrees in case of latlon) 
     472        rMaxY = meta->gds.lat1; // latitude in degrees, to be transformed to meters  
     473        OGRCoordinateTransformation *poTransformLLtoSRS = OGRCreateCoordinateTransformation( &(oLL), &(oSRS) ); 
     474        if ((poTransformLLtoSRS != NULL) && poTransformLLtoSRS->Transform( 1, &rMinX, &rMaxY )) // transform it to meters 
     475        { 
     476            if (meta->gds.scan == GRIB2BIT_2) // Y is minY, GDAL wants maxY 
     477                rMaxY += (meta->gds.Ny - 1) * meta->gds.Dy; // -1 because we GDAL needs the coordinates of the centre of the pixel 
     478        } 
     479        delete poTransformLLtoSRS; 
     480        rPixelSizeX = meta->gds.Dx; 
     481        rPixelSizeY = meta->gds.Dy; 
     482    } 
     483    adfGeoTransform[0] = rMinX; 
     484    adfGeoTransform[3] = rMaxY; 
     485    adfGeoTransform[1] = rPixelSizeX; 
     486    adfGeoTransform[5] = -rPixelSizeY; 
     487 
     488    oSRS.exportToWkt( &(pszProjection) ); 
     489} 
     490 
     491/************************************************************************/ 
     492/*                         GDALRegister_GRIB()                          */ 
    493493/************************************************************************/ 
    494494