WKT for EPSG:4326 with LON 0 to +360 instead of -180 to +180

Question

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:

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

Answer 1

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.

Answer 2

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°.

Answer 3

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

Answer 4

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.

Answer 5

+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