[gdal-dev] question about reading grib2 file

[Matt Funk]

Hi Bill,

thank you. That information helps a lot.
Sorry, as i said, i am pretty new with this stuff.

But it makes sense now.

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
>