[gdal-dev] gdal vs wgrib2 output for grb2 files

[Matt Funk]

Hi Bill,

i was wondering if there had been any more progress on the issue w.r.t
to the discreptancy between wgrib2 and the gdal driver file?
I found the thread which you initiated at:
http://osgeo-org.1803224.n2.nabble.com/gdal-dev-GDAL-Navigation-of-a-Grib2-file-td5740761.html

If you don't mind, can you tell me where that whole thing ended up or
whether it went any further than listed?
Which output did you trust?

Any advice is appreciated...

thanks
matt


On 1/20/2011 3:24 PM, Cassanova, Bill wrote:
> Hi Matt,
>
> Can you be a little more explicit about what you are asking?  Generally speaking GRIB2 files do not contain explicit lat-lon information.
> GRIB2's contain a projection and within that projection the bounding box including the number of points along each axis is contained.
>
> As an experiment run gdalinfo on your grib2 file and you will see output similar to the below:
>
> gdalinfo foo.grib2
> Driver: GRIB/GRIdded Binary (.grb)
> Files: foo.grib2
> Size is 720, 361
> Coordinate System is:
> GEOGCS["Coordinate System imported from GRIB file",
>     DATUM["unknown",
>         SPHEROID["Sphere",6371229,0]],
>     PRIMEM["Greenwich",0],
>     UNIT["degree",0.0174532925199433]]
> Origin = (0.000000000000000,270.000000000000000)
> Pixel Size = (0.500000000000000,-0.500000000000000)
> Corner Coordinates:
> Upper Left  (       0.000,     270.000) (  0d 0'0.01"E,270d 0'0.00"N)
> Lower Left  (   0.0000000,  89.5000000) (  0d 0'0.01"E, 89d30'0.00"N)
> Upper Right (     360.000,     270.000) (360d 0'0.00"E,270d 0'0.00"N)
> Lower Right (     360.000,      89.500) (360d 0'0.00"E, 89d30'0.00"N)
> Center      (     180.000,     179.750) (180d 0'0.00"E,179d45'0.00"N)
> Band 1 Block=720x1 Type=Float64, ColorInterp=Undefined
>   Description = 0[-] SFC="Ground or water surface"
>   Metadata:
>     GRIB_UNIT=[K]
>     GRIB_COMMENT=Temperature [K]
>     GRIB_ELEMENT=TMP
>     GRIB_SHORT_NAME=0-SFC
>     GRIB_REF_TIME=  1295481600 sec UTC
>     GRIB_VALID_TIME=  1296151200 sec UTC
>     GRIB_FORECAST_SECONDS=669600 sec
>     GRIB_PDS_PDTN=0
>     GRIB_PDS_TEMPLATE_NUMBERS=0 0 2 0 96 0 0 0 1 0 0 0 186 1 0 0 0 0 0 255 0 0 0 0 0
>
> If you are trying to read the values contained in Upper Left, Lower Left, etc?
>
> If this is the case you need to get the geo_transform using code similar to the C-code below:
>
> GDALAllRegister();
>
>
>    GDALDatasetH my_data_set;
>
>    int my_x_size, my_y_size;
>    int my_raster_bands;
>
>    GDALRasterBandH my_raster_band;
>
>    double my_geo_transform[6];
>
>    std::string my_projection_ref;
>
>    my_data_set = NULL;
>    my_raster_band = NULL;
>    my_x_size = my_y_size = my_raster_bands = 0;
>
>    ACE_LOG_MSG->log(LM_INFO, "%D %n Opening %s \n", input_grib2_file.c_str());
>
>    my_data_set = GDALOpen(input_grib2_file.c_str(), GA_ReadOnly);
>
>
>    if (NULL == my_data_set)
>    {
>       ACE_LOG_MSG->log(LM_INFO, "%D %n NULL GDALDatasetH pointer for file %s, "
>                        "file invalid, quitting\n", input_grib2_file.c_str());
>       return (1);
>    }
>
>    GDALGetGeoTransform(my_data_set, my_geo_transform);
>
>    my_x_size = GDALGetRasterXSize(my_data_set);
>
>    my_y_size = GDALGetRasterYSize(my_data_set);
>
>    my_raster_bands = GDALGetRasterCount(my_data_set);
>
>    ACE_LOG_MSG->log(LM_INFO,
>                     "%D %n Data Set: x_size = %d y_size = %d raster_bands = %d origin = %f, %f, pixel_size=%f, %f\n",
>                     my_x_size, my_y_size, my_raster_bands, my_geo_transform[0],
>                     my_geo_transform[3], my_geo_transform[1],
>                     my_geo_transform[5]);
>
> >From this information you can calculate the lat-lon of each point within the file.
>
> To answer your other question about wgrib2 -spread....wgrib2 calculates the lat-lon coordinates on the fly and does not store these explicitly in the file...Unless they are put there as a wgrib2 data segment which most people don't do.
>
>
> Hope this helps,
> Bill
>
> -----Original Message-----
> From: gdal-dev-bounces at lists.osgeo.org [mailto:gdal-dev-bounces at lists.osgeo.org] On Behalf Of Matt Funk
> Sent: Thursday, January 20, 2011 4:38 PM
> To: gdal-dev at lists.osgeo.org
> Subject: [gdal-dev] question about reading grib2 file
>
> Hi,
>
> i am somewhat new to gdal and the grib format.
> I need to read in a grib2 file which contains latitude longitude and
> elevation.
> So i do:
> ds = gdal.Open(filename)
> and
> print "ds.GetRasterCount(): %s" % ds.GetRasterCount()
> gives: 1
>
> I can proceed to read in this band, but all the data that i get is the
> elevation.
> Not sure if this is an incredibly stupid question, but how do i access
> the latitude/longitude data?
>
> I know that there is latitude/longitude since when outputting things to
> a csv file via the wgrib2 utility i get lat/lon and elevation?
>
> thanks for any advice
> matt
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>