Deal with 0-360 deg data

[alexander-petkov]

0-360 deg Data

Find a way to deal with data georeferenced in 0-360 deg coordinate space, instead of -180 to 180. Grrr...

Originally posted by @alexander-petkov in https://github.com/alexander-petkov/wfas/issues/3#issuecomment-533283066

[alexander-petkov]

I am creating a new issue, so I can add some notes with my findings in a separate ticket.

Original grid description:

1. Using cdo, the grid from the originally downloaded data is described as follows:

/mnt/cephfs/tmp# cdo griddes gfs.t00z.pgrb2.0p25.f003.grb2
#
# gridID 1
#
gridtype  = lonlat
gridsize  = 40125
xname     = lon
xlongname = longitude
xunits    = degrees_east
yname     = lat
ylongname = latitude
yunits    = degrees_north
xsize     = 321
ysize     = 125
xfirst    = 220
xinc      = 0.25
yfirst    = 22
yinc      = 0.25
cdo griddes: Processed 1 variable ( 0.06s )

2. Using gdalinfo (version 1.11.3, released 2015/09/16):

gdalinfo gfs.t00z.pgrb2.0p25.f003.grb2
Driver: GRIB/GRIdded Binary (.grb)
Files: gfs.t00z.pgrb2.0p25.f003.grb2
Size is 321, 125
Coordinate System is:
GEOGCS["Coordinate System imported from GRIB file",
    DATUM["unknown",
        SPHEROID["Sphere",6371229,0]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]
Origin = (219.875000000000000,53.125000000000000)
Pixel Size = (0.250000000000000,-0.250000000000000)
Corner Coordinates:
Upper Left  (     219.875,      53.125) (219d52'30.00"E, 53d 7'30.00"N)
Lower Left  (     219.875,      21.875) (219d52'30.00"E, 21d52'30.00"N)
Upper Right (     300.125,      53.125) (300d 7'30.00"E, 53d 7'30.00"N)
Lower Right (     300.125,      21.875) (300d 7'30.00"E, 21d52'30.00"N)
Center      (     260.000,      37.500) (260d 0' 0.00"E, 37d30' 0.00"N)

3. Using gdalinfo (version 2.4.0, released 2018/12/14): Interestingly enough, this much newer version reports lat/lon extents in -180/180 and -90/90 deg space:

gdalinfo gfs.t00z.pgrb2.0p25.f003 
Driver: GRIB/GRIdded Binary (.grb, .grb2)
Files: gfs.t00z.pgrb2.0p25.f003
Size is 321, 125
Coordinate System is:
GEOGCS["Coordinate System imported from GRIB file",
    DATUM["unknown",
        SPHEROID["Sphere",6371229,0]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]
Origin = (-140.125000000000000,53.125000000000000)
Pixel Size = (0.250000000000000,-0.250000000000000)
Corner Coordinates:
Upper Left  (-140.1250000,  53.1250000) (140d 7'30.00"W, 53d 7'30.00"N)
Lower Left  (-140.1250000,  21.8750000) (140d 7'30.00"W, 21d52'30.00"N)
Upper Right ( -59.8750000,  53.1250000) ( 59d52'30.00"W, 53d 7'30.00"N)
Lower Right ( -59.8750000,  21.8750000) ( 59d52'30.00"W, 21d52'30.00"N)
Center      (-100.0000000,  37.5000000) (100d 0' 0.00"W, 37d30' 0.00"N)

[alexander-petkov]

Re-writing the lat/lon extents for the data

gdalwarp can be used to rewrite the data extent around a 0 deg prime meridian: gdalwarp -t_srs WGS84 gfs.t00z.pgrb2.0p25.f069.grb2 /tmp/output180.tif --config CENTER_LONG 0 -wo SOURCE_EXTRA=100 -overwrite

My problems with using gdal are as follows.

1. I don't have gdal on the VM which pulls the data.
2. If gdal is used, the time stamp might have to be encoded in the file name. Actually, GDAL 2.4 supports writing to Grib format.
3. I already compiled cdo, which seems to be pretty capable with NetCDF and Grib data, and I already use it during data download to derive other datasets.

This forum discussion has been very informative. I tried to alter the extent using cdo sellonlatbox , which did not work. I tried the next suggesstion, which is to overwrite the grid definition:

you can do a 360 degree shift of the longitudes by overwriting the grid description with setgrid if you have a regular lonlat grid which crosses neither the prime meridian nor the date line:

cdo griddes ifile > mygrid
edit mygrid and set xfirst to the new value
cdo setgrid,mygrid ifile ofile

Changing the xfirst value to -140 should solve my issue, right? Wrong, the grid in the new file remained unaltered.

I started using different values, to see if rewriting the grid had any effect. xfirst = -360 put the left boundary at 0 degrees? Using even lower values yelded undesired outcome, with bouding box to the effect of 2000 degrees longitude.

[alexander-petkov]

What worked

These two settings in the grid description allowed me to correct the problem:

1. xunits = degrees_west instead of degrees_east
2. xfirst = -500 ??? Now this is puzzling, because a value of -140 should be correct. Do I need to wrap the data around the globe a time and a half? Very strange:

cdo sinfon gfs.t00z.pgrb2.0p25.f003.grb2.new 
   File format : GRIB2
    -1 : Institut Source   Ttype    Levels Num    Points Num Dtype : Parameter name
     1 : NCEP     unknown  instant       1   1     40125   1  P9   : 2t            
   Grid coordinates :
     1 : lonlat                   : points=40125 (321x125)
                              lon : -140 to -60 by 0.25 degrees_east
                              lat : 22 to 53 by 0.25 degrees_north

gdalinfo shows that Earth got a whole lot bigger:

gdalinfo gfs.t00z.pgrb2.0p25.f003.grb2.new 
Driver: GRIB/GRIdded Binary (.grb)
Files: gfs.t00z.pgrb2.0p25.f003.grb2.new
Size is 321, 125
Coordinate System is:
GEOGCS["Coordinate System imported from GRIB file",
    DATUM["unknown",
        SPHEROID["Sphere",6367470,0]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]
Origin = (-2007.608647999999903,53.125000000000000)
Pixel Size = (0.250000000000000,-0.250000000000000)
Corner Coordinates:
Upper Left  (   -2007.609,      53.125) (Invalid angle, 53d 7'30.00"N)
Lower Left  (   -2007.609,      21.875) (Invalid angle, 21d52'30.00"N)
Upper Right (   -1927.359,      53.125) (Invalid angle, 53d 7'30.00"N)
Lower Right (   -1927.359,      21.875) (Invalid angle, 21d52'30.00"N)
Center      (   -1967.484,      37.500) (Invalid angle, 37d30' 0.00"N)

Configuring this new file as a coverage in Geoserver:

New file displayed as a WCS layer in QGIS:

[alexander-petkov]

Self-compiled cdo did not work the same on manager1, where data collection is performed via scheduled cron tasks.

I had compiled cdo on a VM that is running Ubuntu 16.x release. Rewriting the grid on it worked fine., but not so on manager1, which has the same ceph pool mounted.

I solved this by installing cdo via apt. It took additional 1.3G of disk space. A bit excessive.

Also, use a newer gdal: version. 2.4.0 on Ubuntu 19 seems to work as expected. Degrees are still converted to Celsius though...

[alexander-petkov]

The latest grid definition that works with cdo:

#
# gridID 1
#
gridtype  = latlon
gridsize  = 40125
xname     = lon
xlongname = "longitude"
xunits    = "degrees_west"
yname     = lat
ylongname = "latitude"
yunits    = "degrees_north"
xsize     =  321
ysize     =  125
xfirst    =  -140
xinc      =  -0.25
yfirst    =  22
yinc      =  0.25
scanningMode = 64

Update: this works with cdo, and also an up to date GDAL version. However, Geoserver goes williwompus when computing the data extent:

The struggle continues...

[alexander-petkov]

As of now, I have solved this by converting downloaded GFS data to nc format.