5

I have GRIB2 files that seem to be in a projection similar to EPSG: 4326, except that instead of longitude between -180 and +180 it is between 0 and 360.

My intention is to get these files into geoserver and allow the data to be queried and reprojected at will through the WMS. I can get these files into geoserver with the GRIB extension, but I get a null pointer error when I attempt to reproject to mercator. I assume that this error occurs because data lies outside the -180 to +180 longitude bound.

I had attempted to just reproject the GRIB files with gdal_translate, using the PROJ4 +lon_wrap=180option, but that operation takes forever and produces a far larger output file for some reason.

From what I understand, I can fix this situation by creating a custom projection in user_projections/epsg.propertiesthat describes the CRS of these GRIB files. A custom projection for geoserver must be in the WKT format.

TL;DR: How do I write a WKT that describes a lat/lon coordinate system similar to EPSG: 4326 but with longitude ranging from 0 to +360 instead of -180 to +180?

Here is the gdalinfo for one of these GRIB files. The file depicts the continental United States:

Size is 2847, 1517
Coordinate System is:
GEOGCS["Coordinate System imported from GRIB file",
    DATUM["unknown",
        SPHEROID["Sphere",6371229,0]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]
Origin = (229.982699929725896,55.001944129287601)
Pixel Size = (0.024600140548138,-0.023088258575198)
Corner Coordinates:
Upper Left  (     229.983,      55.002) (229d58'57.72"E, 55d 0' 7.00"N)
Lower Left  (     229.983,      19.977) (229d58'57.72"E, 19d58'37.40"N)
Upper Right (     300.019,      55.002) (300d 1' 9.48"E, 55d 0' 7.00"N)
Lower Right (     300.019,      19.977) (300d 1' 9.48"E, 19d58'37.40"N)
Center      (     265.001,      37.490) (265d 0' 3.60"E, 37d29'22.20"N)

gdalsrsinfo:

PROJ.4 : '+proj=longlat +a=6371229 +b=6371229 +no_defs '

OGC WKT :
GEOGCS["Coordinate System imported from GRIB file",
    DATUM["unknown",
        SPHEROID["Sphere",6371229,0]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]

Screenshot from QGIS: QGIS screenshot

Edit: This is not a duplicate. I have seen the answer from the referenced question, and it is not satisfactory. The operation takes too long to complete and produces a file that is entirely too large.

Edit #2: Added a link to the GRIB file: http://furlender.com/ECMWF_2017052412.grib2

  • There's no parameter that I know of in the WKT1 or 2 specifications. Could you re-georeference the file? – mkennedy Sep 21 '17 at 19:12
  • Really... thats too bad. You can define all these complex projections but you can't just offset the latitude by 180? – Mike Furlender Sep 21 '17 at 19:56
  • Can you share a download link for the file? – dbaston Sep 22 '17 at 16:35
  • furlender.com/ECMWF_2017052412.grib2 @dbaston – Mike Furlender Sep 22 '17 at 17:14
  • 1
    @rensa Nope - I ended up just warping it with gdalwarp. From what I can tell, there's no way to produce a WKT that describes this. – Mike Furlender Feb 21 '18 at 21:21
4

You can simply build a vrt file around your source grib with

gdal_translate -of VRT -a_ullr -130.0173001 55.0019441 -59.9870389 19.9746766 ECMWF_2017052412.grib2 ECMWF180.vrt

which is just 47kB large. The extent is calculated by hand. Similar to How to reproject raster from 0 360 to -180 180 with cutting 180 meridian you might as well run

gdal_translate -of VRT -b 1 ECMWF_2017052412.grib2 temp.vrt
gdalwarp -overwrite -t_srs WGS84 temp.vrt ECMWF180.tif  --config CENTER_LONG 0
gdalinfo ECMWF180.tif

to extract it from the first band. The grib file has 41 bands with a fairly high compression; converting all to tif uses over 1GB and 2 min calculation.

You might add -a_srs EPSG:4326 to the simple solution if you don't mind the difference between the source sphere and WGS84 ellipsoid. The simple solution will work unless you hit the 0 or 180° meridians, while the second will work always.

  • This seems like a workable solution. Thanks a ton! – Mike Furlender Sep 23 '17 at 19:27
2

You might want to give wgrib2 a try. From this file of wgrib2 tips & tricks is the following:

(42) I want to convert a grid to -180 to 180 grid.

Assuming that you have IPOLATES installed and IPOLATES supports the grid, then you can do:

wgrib2 In.grb -new_grid_winds earth -new_grid latlon 180:361:1 -90:181:1 out.grb

If you don't want an 1x1 degree grid, you would have to change the grid definition.

In my experience, "assuming that you have IPOLATES installed" is equivalent to "assuming that you compiled wgrib2 yourself," but perhaps there are binaries with IPOLATES out there somewhere.

  • This looks like a potential solution to my problem, but when I use a 1x1 grid I get a tiny file with a very pixelated grid. How do I find out the dimensions of my original grid? I thought it might be "Size is 2847, 1517" from gdalinfo, but using 180:361:1517 -90:181:2847 out.grb gave me errors. – Mike Furlender Sep 21 '17 at 18:38
  • I think the example from the FAQ might have a typo. The syntax is min:ncells:resolution. Sounds like your inputs are global at 1/8-degree resolution, so you could do -new_grid latlon -180:2880:0.125 -90:1440:0.125. – dbaston Sep 21 '17 at 18:55
  • That gets me closer, but when I do gdalinfo on the resulting file I get Upper Left ( 179.938, 99.562) (179d56'15.00"E, 99d33'45.00"N) Lower Left ( 179.9375000, -90.0625000) (179d56'15.00"E, 90d 3'45.00"S) Upper Right ( 535.812, 99.562) (Invalid angle, 99d33'45.00"N) Lower Right ( 535.812, -90.062) (Invalid angle, 90d 3'45.00"S) Center ( 357.875, 4.750) (357d52'30.00"E, 4d45' 0.00"N) – Mike Furlender Sep 21 '17 at 19:01
  • This page has a better description of the latlon parameter. You need start point (which is not 180,-90), number of cells, cell size so it's basically regeoreferencing the data if you use the same cell size. – mkennedy Sep 21 '17 at 20:30
  • this was my latest attempt with wgrib2: wgrib2 ECMWF_2017052412.grib2 -new_grid_winds earth -new_grid latlon -130.0173001:2847:0.0174532925199433 19.9770559:1517:0.125 out.grb – Mike Furlender Sep 22 '17 at 19:49
1

I just ran into this with similar GRIB data... both MRMS radar data and AirNow air quality data are projected in 0-360° longitude, I'm assuming to save space by using unsigned integers.

Using the "+pm" parameter of PROJ fixed this:

+proj=longlat +ellps=WGS84 +pm=-360 +datum=WGS84 +no_defs

For a MapServer layer, this worked out to:

PROJECTION
###    "init=epsg:4326"
       "+proj=longlat +ellps=WGS84 +pm=-360 +datum=WGS84 +no_defs"
END

I've not used GeoServer, but this might work for the WKT description:

GEOGCS["4326 plus 360",
    DATUM["unknown",
        SPHEROID["Sphere",6371229,0]],
    PRIMEM["Custom",-360],
    UNIT["degree",0.0174532925199433]]

The reason for the "-360" is that the coordinates listed in the original question (as well as my own) have an Upper Left of ~230°, when that should actually be translated/rotated to -130°. i.e.: (230-360) = -130°.

1

Although these solutions might have worked for some, in my case, I needed to do rendering transformations of .grib2 rasters on the fly. Geoserver WMS reprojection worked well for the raster itself but the created contours and wind arrows were created only for 0 - 360 longitude and not for -180,0.

Although I wanted to avoid using gdalwarp in my workflow, I had to change my mind and had to warp the grib2 to .tif, which finally works.

gdalwarp -t_srs WGS84 input360.grib2 output180.tif --config CENTER_LONG 0 -wo SOURCE_EXTRA=100 -overwrite

0

+pm=-360 aligned the data for me in mapserver but still wouldn't render properly. I had to replace it with +lon_wrap=180 to get it to work

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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