grib2 to raster? (grib2 -> raster2pgsql -> PostGIS 2.0 -> GeoServer)

Question

I'm looking to implement import of US National Weather Service data for doing various sorts of models (www.soils.wisc.edu/uwex_agwx if you're curious about what we're doing with the data).

The weather data come in Grib2 files, and I'm fiddling about with PostGIS 2.0 in an attempt to make them accessible to both my own modeling software and to put up visualizations via GeoServer.

I can get it to work with some files, or at least I get SQL out that seems to run OK, creates a table and loads data into it. But I cannot visualize the resulting table in GeoServer, so I know that something is wacky along the way.

So my questions -- pardon their newbie-tude:

1) Anybody done anything like this here before? I've seen folks trying to build things out of GRiBs, but not rasters.

2) I have been using the Rapid Update Cycle forecast data, in which coordinates for the continental US look something like e.g. X=238.1234, Y=10.1234. From this I'm assuming that I can use EPSG:4326 as an SRID -- am I correct or crazy here? Anyone know how I might induce the correct SRID code to use from the coordinates?

If anyone would like to play, here's a 35MB grib2 file. I used this command line to create a table and insert those data into it: raster2pgsql -s 4326 -C -M ruc2.t00z.bgrb13anl.grib2 grib_tbl

NOTE: I cross-posted this to StackOverflow, and got the suggestion there of running gdalinfo on the GRIB2 file. More as it develops. I will ensure that if this gets answered so as to close out the question in either forum, I WILL SO POST IN BOTH PLACES!

Answer 1

Looking at the data file you have posted, I can see right away that your assumption is incorrect.

You can't just assume you have the right coordinate system, you have to know what your system is. You can get this from running gdalinfo, as you mentioned. That gives us:

Driver: GRIB/GRIdded Binary (.grb)
Files: /tmp/ruc2.t00z.bgrb13anl.grib2
Size is 451, 337
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["Coordinate System imported from GRIB file",
        DATUM["unknown",
            SPHEROID["Sphere",6371229,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",25],
    PARAMETER["standard_parallel_2",25],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",265],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0]]
Origin = (-3332155.288903323933482,6830293.833488883450627)
Pixel Size = (13545.000000000000000,-13545.000000000000000)
Corner Coordinates:
Upper Left  (-3332155.289, 6830293.833) (139d51'22.04"W, 54d10'20.71"N)
Lower Left  (-3332155.289, 2265628.833) (126d 6'34.06"W, 16d 9'49.48"N)
Upper Right ( 2776639.711, 6830293.833) ( 57d12'21.76"W, 55d27'10.73"N)
Lower Right ( 2776639.711, 2265628.833) ( 68d56'16.73"W, 17d11'55.33"N)
Center      ( -277757.789, 4547961.333) ( 98d 8'30.73"W, 39d54'5.40"N)

What this tells me is that we don't have a common coordinate system. (How do I know? The unknown datum, the spheroid as "Sphere" are the big clues.) That means it's not going to be able to be determined by the raster2pgsql command (it will try, but if it fails, assumes no coordinate system.)

But, GDAL is capable of dealing with this type of bespoke coordinate system, so the way to handle it is to have GDAL transform it to another system that PostGIS understands, and import that new file:

gdalwarp -t_srs EPSG:4326 ruc2.t00z.bgrb13anl.grib2 transformed.grib2

That's simply asking GDAL to write a new file, but change the projection to WGS84. Now you should be able to import that into postgres using the command you have above.

FWIW, that data file has 750 bands, which seems to eat a lot of memory on my machine on import. You may want to limit the import to only the bands you want. Here's how imported only band 1.

raster2pgsql -I -b 1 -C -M /tmp/transformed.grib2  grib_test |  psql raster_test

Good luck.