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=180
option, 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.properties
that 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
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.
No comments:
Post a Comment