GRIB-2/ICON geolocation unknown or invalid for western hemisphere

[gerritholl]

Describe the bug

When reading ICON data in GRIB-2 format using the grib reader, and showing it along with coastlines, there are no coastlines for the western hemisphere. Geolocation appears to accept only longitudes between 0 and 180 degrees.

To Reproduce

from satpy import Scene
from glob import glob
sc = Scene(filenames={"grib":
    glob("/media/nas/x21308/scratch/ICON/*12:00:00Z*000.grib")})
sc.load(["t"])
sc.show("t", overlay={"coast_dir": "/media/nas/x21308/shp/", "color": "red"})

Expected behavior

I expect a world map with the field displayed and worldwide coastlines. This map would probably be centered on the meridian.

Actual results

The map is centered on the antimeridian and coastlines are shown for the eastern hemisphere only.

Environment Info:

• OS: openSUSE 15.0
• Satpy Version: 0.21.1.dev184+g836c657d
• PyResample Version: 1.15.0+63.g3b584e6
• Readers and writers dependencies (when relevant):

Readers and writers dependencies

Readers                                                                                                                                                                                                                                                                          [24/1823]
=======
abi_l1b:  ok
abi_l1b_scmi:  ok
abi_l2_nc:  ok
acspo:  ok
agri_l1:  ok
ahi_hrit:  ok
ahi_hsd:  ok
ami_l1b:  ok
amsr2_l1b:  ok
avhrr_l1b_aapp:  ok                                                                                                                                                                                                                                                                      
avhrr_l1b_eps:  ok
avhrr_l1b_gaclac:  cannot find module 'satpy.readers.avhrr_l1b_gaclac' (No module named 'pygac')
avhrr_l1b_hrpt:  cannot find module 'satpy.readers.hrpt' (No module named 'pygac')
caliop_l2_cloud:  ok
clavrx:  ok
cmsaf-claas2_l2_nc:  ok
electrol_hrit:  ok
fci_l1c_fdhsi:  ok
generic_image:  ok
geocat:  ok
ghrsst_l3c_sst:  ok
glm_l2:  ok
goes-imager_hrit:  ok
goes-imager_nc:  ok
grib:  ok
hsaf_grib:  ok
hy2_scat_l2b_h5:  ok
iasi_l2:  ok
iasi_l2_so2_bufr:  ok
jami_hrit:  ok
li_l2:  cannot find module 'satpy.readers.li_l2' (No module named 'h5netcdf')
maia:  ok
mersi2_l1b:  ok
mimicTPW2_comp:  ok
modis_l1b:  ok
modis_l2:  ok
msi_safe:  cannot find module 'satpy.readers.msi_safe' (No module named 'glymur')
mtsat2-imager_hrit:  ok
nucaps:  ok
nwcsaf-geo:  ok
nwcsaf-msg2013-hdf5:  ok
nwcsaf-pps_nc:  ok
olci_l1b:  ok
olci_l2:  ok
omps_edr:  ok
safe_sar_l2_ocn:  ok
sar-c_safe:  ok
scatsat1_l2b:  ok
seviri_l1b_hrit:  ok
seviri_l1b_icare:  ok
seviri_l1b_native:  ok
seviri_l1b_nc:  ok
seviri_l2_bufr:  ok
seviri_l2_grib:  ok
slstr_l1b:  ok
slstr_l2:  ok
tropomi_l2:  ok
vaisala_gld360:  ok
viirs_compact:  ok
viirs_edr_active_fires:  ok
viirs_edr_flood:  ok                                                             
viirs_l1b:  ok         
viirs_sdr:  ok
virr_l1b:  ok  
                        
Writers           
=======      
/data/gholl/miniconda3/envs/py38/lib/python3.8/site-packages/pyninjotiff/tifffile.py:154: UserWarning: failed to import the optional _tifffile C extension module.
Loading of some compressed images will be slow.
Tifffile.c can be obtained at http://www.lfd.uci.edu/~gohlke/
  warnings.warn(
cf:  ok          
geotiff:  ok        
mitiff:  ok          
ninjotiff:  ok        
scmi:  ok         
simple_image:  ok  
                   
Extras        
======       
cartopy:  ok
geoviews:  ok

Additional context

The AreaDefinition is:

Area ID: on-the-fly grib area
Description: on-the-fly grib area
Projection ID: on-the-fly grib area
Projection: {'R': '6371229', 'lat_0': '0', 'lat_ts': '0', 'lon_0': '0', 'no_defs': 'None', 'proj': 'eqc', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}                                                                                                                           
Number of columns: 1440
Number of rows: 721
Area extent: (-13899.8654, -10021802.9758, 40017712.576, 10021802.9758)

The get_xy_from_lonlat method appears to accept only longitudes between 0 and 180. For any longitude smaller than 0 or larger than 180, it throws an exception ValueError: Point outside area (or issues NaN when passed an array).

[djhoese]

Looks like the data is probably 0 to 360 longitude instead of -180 to 180. There is some special handling this for lon/lat projections, but not eqc. There is even some +over PROJ parameter usage for the cyl projection that I was seeing in some of the NCEP/NOAA files and is equivalent to eqc. The code could probably be updated to handle cyl and eqc as the same thing and this would "just work". Maybe.

[gerritholl]

What is a "lon/lat" projection? EQC = equirectangular = Plate Carrée?

The projparams in this GRIB file are {'a': 6371229, 'b': 6371229, 'proj': 'cyl'}.

[gerritholl]

I confirm that the longitudes go from 0 to 360, at least that's what msg.latlons() gives.

[djhoese]

"Lon/lat" is how I refer to +proj=longlat (technically PROJ understands lonlat, latlong, and maybe latlon). This is measured in degrees (or radians). The eqc projection is measured in meters.

[gerritholl]

When you say "special handling code", do you mean this part?

satpy/satpy/readers/grib.py

Lines 145 to 147 in ab17a71

	for lon_param in ['lon_0', 'lon_1', 'lon_2']:
	if proj_params.get(lon_param, 0) > 180:
	proj_params[lon_param] -= 360

That seems to be called in any case. I don't quite understand what happens in the rest of this method.

[djhoese]

Ha yes, sorry, I meant to link directly to that but forgot the link. There's this part:

satpy/satpy/readers/grib.py

Line 167 in ab17a71

proj_params['over'] = True

and this part:

satpy/satpy/readers/grib.py

Line 185 in ab17a71

lons[lons > 180] -= 360

All to handle data not being in the right order or at least not in the expected order. Since everyone uses GRIB the way they want, it was really hard to come up with something simple that worked for every file. This was as far as I got for the files I needed to support. I think if the below lines were updated:

satpy/satpy/readers/grib.py

Lines 149 to 150 in ab17a71

	if proj_params['proj'] == 'cyl':
	proj_params['proj'] = 'eqc'

to:

        if proj_params['proj'] == 'cyl':
            proj_params['proj'] = 'eqc'
        if proj_params['proj'] == 'eqc':

Your case might work.

[gerritholl]

Thanks for the suggestion. Unfortunately the image still have half of the coastlines missing just like in the image I posted before.

[djhoese]

That could be a separate pycoast issue that I think there is an issue for but can't find it right now. Depends what the AreaDefinition is now.

[gerritholl]

The AreaDefinition is the same as before:

Area ID: on-the-fly grib area
Description: on-the-fly grib area
Projection ID: on-the-fly grib area
Projection: {'R': '6371229', 'lat_0': '0', 'lat_ts': '0', 'lon_0': '0', 'no_defs': 'None', 'proj': 'eqc', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
Number of columns: 1440
Number of rows: 721
Area extent: (-13899.8654, -10021802.9758, 40017712.576, 10021802.9758)

and it still won't let me get the xy for longitudes -1 or 181.

[gerritholl]

It does enter the block starting with the comment # wrap around.

[gerritholl]

Could it be a problem in pyresample? I found this surprising:

ar1 = AreaDefinition("test", "test", "test", {'a': 6371229, 'b': 6371229, 'proj': 'eqc', 'over': True}, 1440, 721, (-13899.865431068225, -10021802.97580019, 40017712.576045424, 10021802.97580019))
ar2 = AreaDefinition("test", "test", "test", {'a': 6371229, 'b': 6371229, 'proj': 'eqc', 'over': False}, 1440, 721, (-13899.865431068225, -10021802.97580019, 40017712.576045424, 10021802.97580019))
print(ar1 == ar2)

which gives True.

It looks like the over parameter has no effect?

[djhoese]

That's likely PROJ removing it. Older versions of PROJ used it as a hack. The best long term solution may be to cut the array up and reorder it to correspond to the projection. Right now the data is going from 0 to 360. The projection space though goes from -180 to 180. So in the AreaDefinition you have a non-contiguous Area: 0 to 180, -180 to 0. In meters this is probably something like ~0 to 20000000, -20000000 to 0. If the two halves of the data switched places in the array then it would work. I'm not sure there is a standard way in PROJ nowadays to be in a 0 to 360 degree space.

Additionally, pyresample may have some limitations right now regarding lon/lat -> lon/lat (see pytroll/pyresample#267 and point 1).

[gerritholl]

So probably this affect any GRIB data that has longitudes 0 to 360, not only ICON?

[djhoese]

I don't know what ICON is, but yes probably. Now that I think about it this probably effects the GRIB files that I was using for my SIFT project with newer PROJ versions.

[gerritholl]

ICON is the global model developed and used by the German Weather Service (see https://www.dwd.de/SharedDocs/downloads/DE/modelldokumentationen/nwv/icon/icon_dbbeschr_aktuell.pdf) but I suspect the problem described here shouldn't be specific to ICON (i.e. there's probably nothing wrong with the ICON files; the NWCSAF software can interpret them apparently correctly).

[simonrp84]

I have noticed a similar problem for ECMWF forecasts, which also use 0-360 for longitude.

[djhoese]

I've run in to this issue affecting my SIFT application. I'll see if I can come up with a fix that works with newer versions of PROJ/pyproj.