2

I'm working on a subsetted grib file that I retrieved using this URL: http://nomads.ncep.noaa.gov/cgi-bin/filter_gfs_0p25.pl?file=gfs.t12z.pgrb2.0p25.f001&lev_925_mb=on&var_RH=on&subregion=&leftlon=-126&rightlon=-66&toplat=50&bottomlat=24&dir=%2Fgfs.2017020912

First I opened the grib file in zyGrib to confirm that the file has the right data, i.e., relative humidity is appearing for lower 48 US states at an altitude of 925hPa. So far so good.

Next, I try loading the grib file into postgres:

raster2pgsql -s 4326 gfs.t12z.pgrb2.0p25.f001 | psql

Then, I confirm the new raster exists, has one band, and SRID=4326:

SELECT ST_Height(rast) As height,
ST_Width(rast) As width,
ST_SRID(rast) As srid,
ST_NumBands(rast) As num_bands,
ST_BandPixelType(rast,1) As btype
FROM public."gfs.t12z.pgrb2.0p25";

Next, I try querying a point located in the USA:

SELECT ST_Value(rast, ST_SetSRID(ST_Point(-87.9080059, 41.9773201), 4326))
FROM public."gfs.t12z.pgrb2.0p25";

I'm getting the following warning:

NOTICE: Attempting to get pixel value with out of range raster coordinates: (-1288, 32)

And a blank value is returned in the resultset.

So next, I try gdalwarp as suggested in this answer: https://gis.stackexchange.com/a/33524/12291

gdalwarp -t_srs EPSG:4326 gfs.t12z.pgrb2.0p25.f001 warped.grib2

Then import:

raster2pgsql -s 4326 warped.grib2 | psql

Same as before, confirm new raster exists, has one band, and SRID=4326. Results match last attempt.

SELECT ST_Height(rast) As height,
ST_Width(rast) As width,
ST_SRID(rast) As srid,
ST_NumBands(rast) As num_bands,
ST_BandPixelType(rast,1) As btype
FROM public.warped;

Querying a point located in the USA from the new warped raster:

SELECT ST_Value(rast, ST_SetSRID(ST_Point(-87.9080059, 41.9773201), 4326))
FROM public.warped;

Gives the same result ... warning message and blank value in resultset :-(

NOTICE: Attempting to get pixel value with out of range raster coordinates: (-1288, 32)

What am I doing wrong here?

EDIT: here's the gdalinfo for this grib file:

gdalinfo gfs.t12z.pgrb2.0p25.f001
Driver: GRIB/GRIdded Binary (.grb)
Files: gfs.t12z.pgrb2.0p25.f001
Size is 241, 105
Coordinate System is:
GEOGCS["Coordinate System imported from GRIB file",
    DATUM["unknown",
        SPHEROID["Sphere",6371229,0]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]
Origin = (233.875000000000000,50.125000000000000)
Pixel Size = (0.250000000000000,-0.250000000000000)
Corner Coordinates:
Upper Left  (     233.875,      50.125) (233d52'30.00"E, 50d 7'30.00"N)
Lower Left  (     233.875,      23.875) (233d52'30.00"E, 23d52'30.00"N)
Upper Right (     294.125,      50.125) (294d 7'30.00"E, 50d 7'30.00"N)
Lower Right (     294.125,      23.875) (294d 7'30.00"E, 23d52'30.00"N)
Center      (     264.000,      37.000) (264d 0' 0.00"E, 37d 0' 0.00"N)
Band 1 Block=241x1 Type=Float64, ColorInterp=Undefined
  Description = 92500[Pa] ISBL="Isobaric surface"
  Metadata:
    GRIB_COMMENT=Relative Humidity [%]
    GRIB_ELEMENT=RH
    GRIB_FORECAST_SECONDS=3600 sec
    GRIB_PDS_PDTN=0
    GRIB_PDS_TEMPLATE_NUMBERS=1 1 2 0 96 0 0 0 1 0 0 0 1 100 0 0 1 105 84 255 0 0 0 0 0
    GRIB_REF_TIME=  1486641600 sec UTC
    GRIB_SHORT_NAME=92500-ISBL
    GRIB_UNIT=[%]
    GRIB_VALID_TIME=  1486645200 sec UTC

EDIT #2: Here is the information zyGrib shows for this file. Somehow it is generating the correct area for this file despite missing datum in the source file.

File : gfs.t12z.pgrb2.0p25-7.f001
Size : 31811 bytes

Weather center 7 - Model 96 - Grid 0
1 dates:
    from Thu 2017-02-09 13:00 UTC
    to Thu 2017-02-09 13:00 UTC

Avalaible data :
* Relative humidity: 925hpa

Grid : 241x105=25305 points : 0.25°x0.25°
Area :24.00°N 126.00°W  ->  50.00°N 066.00°W

Reference date: Thu 2017-02-09 13:00 UTC 
6
  • using gdalinfo or qgis, origin shows up as 233.875, 50.125 which isn't valid in wgs84, which might explain shift in x coordinates. (Shows up in QGIS at right latitude, but shifted off the right of the antimeridian). I don't use GRIB files so I'm not sure how you might correct that.
    – Steven Kay
    Commented Feb 9, 2017 at 20:59
  • 1
    Thanks @StevenKay ... that reminded me I forgot to include gdalinfo, have added that to the question.
    – boozedog
    Commented Feb 9, 2017 at 21:51
  • Yeh you have to fix the original dataset, the coordinates of your dataset range between 233.875 to 294.125 latitudes and your trying to query a point that has a lat of -87.9
    – Liam G
    Commented Feb 9, 2017 at 22:50
  • Thanks @Freight_Train ... any suggestion on how I'd fix it? Also, I wonder how come zyGrib seems to handle this file just fine? It shows the raster in exactly the right spot.
    – boozedog
    Commented Feb 9, 2017 at 23:17
  • Do you know what datum zyGrib is using to display the data? I notice that it is unknown in gdal_info
    – Liam G
    Commented Feb 9, 2017 at 23:19

2 Answers 2

2

Got this figured out.

Squinting at the coordinates a bit, I realized that the longitude in this grib file is expressed as 0 to 360 degrees with Greenwich as the prime meridian.

But latitude is expressed as usual, -90 to 90.

Anyway the solution is to subtract the degrees west longitude from 360.

So, the following query:

SELECT ST_Value(rast, ST_SetSRID(ST_Point(-87.9080059, 41.9773201), 4326)) 
FROM public."gfs.t12z.pgrb2.0p25";

Gets rewritten as follows:

/* 360 - 87.9080059 = 272.0919941 */
SELECT ST_Value(rast, ST_SetSRID(ST_Point(272.0919941, 41.9773201), 4326)) 
FROM public."gfs.t12z.pgrb2.0p25";

Also note that I did NOT need to use gdalwarp. Importing the grib as-is using raster2pgsql -s 4326 works fine.

1
  • BTW thanks to the commenters above for providing the hints I needed to figure this out :-)
    – boozedog
    Commented Feb 11, 2017 at 20:54
1

You can also use ST_Shift_Longitude postgis function instead of calculate it.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.