Grib2 to PostGIS raster -- anyone get this to work?

3.6k views Asked by At

I have an application for which I need to import U.S. National Weather Service surface analyses, which are distributed as grib2 files. I want to pull those into PostGIS 2.0 rasters, do some calculations and modeling, and display the data and model results in GeoServer.

Since grib2 is a GDAL-supported format, the supplied raster2pgsql utility should be able to slurp a grib2 right into PostGIS-compatible SQL, and once it's there, GeoServer ought to be able to handle it. However, I'm running into problems which have no obvious solutions -- not obvious to me, at any rate! Raster2pgsql runs, apparently without errors, producing SQL, and running the SQL creates what looks very much like a raster. But GeoServer can't display it -- the bounds, in particular, come out looking weird (0,0 -1,-1) and "preview layer" just throws a NullPointerException.

Has anyone been down this road already? I've got issues as basic as not knowing what the SRID should be for the data (4326, perhaps?). I don't expect anyone to debug my problems for me but if someone has already got this toolchain working, or part of it, I can plug known-good things in and see what I can discover.

TIA,

rw

Updated: Per Mike, here is the coordinate-system stuff from one of the files; I elided the other 749 bands in the output from "gdalinfo". Note that the filename is different -- I found out by running "gdalinfo" on my original file that something was wrong with it, gdalinfo couldn't read it. New (35MB!) file here.

Gdalinfo output:

Driver: GRIB/GRIdded Binary (.grb)
Files: 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)
Band 1 Block=451x1 Type=Float64, ColorInterp=Undefined
  Description = 1[-] HYBL="Hybrid level"
  Metadata:
    GRIB_UNIT=[Pa]
    GRIB_COMMENT=Pressure [Pa]
    GRIB_ELEMENT=PRES
[Etc., Etc., for all 750 bands]
2

There are 2 answers

0
Rick Wayne On BEST ANSWER

Got an excellent answer to my problem here. Putting it in as a separate answer.

He recommended using gdalwarp to pull the GRIB2 file into a known SRID, thus:

gdalwarp -t_srs EPSG:4326 original_file.grib2 4326_file.grib2

Then, raster2pgsql works just fine, e.g.

raster2pgsql -M -a 4326_file.grib2 some_sql.sql
0
Mikel G. Gainza On

I hope this helps, at least those comming to this thread.

Bear in mind that GeoServer, while being capable of loading Raster data from PostGIS, the default PostGIS "importing" module is ONLY available for vector data, that's why you get those odd bounds (-1 -1 0 0).

You'll have to add ImageMosaicJDBC plugin to your geoserver installation, follow steps here!

http://docs.geoserver.org/latest/en/user/tutorials/imagemosaic-jdbc/imagemosaic-jdbc_tutorial.html