Wednesday, 3 July 2019

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



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



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

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...