weixin_39698217
weixin_39698217
2021-01-01 02:00

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

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

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

tmp9v5spg4z

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

该提问来源于开源项目:pytroll/satpy

  • 点赞
  • 写回答
  • 关注问题
  • 收藏
  • 复制链接分享
  • 邀请回答

17条回答

  • weixin_39698217 weixin_39698217 4月前

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

    点赞 评论 复制链接分享
  • weixin_39732027 weixin_39732027 4月前

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

    点赞 评论 复制链接分享
  • weixin_39963255 weixin_39963255 4月前

    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.

    点赞 评论 复制链接分享
  • weixin_39963255 weixin_39963255 4月前

    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.

    点赞 评论 复制链接分享
  • weixin_39698217 weixin_39698217 4月前

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

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

    点赞 评论 复制链接分享
  • weixin_39698217 weixin_39698217 4月前

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

    点赞 评论 复制链接分享
  • weixin_39963255 weixin_39963255 4月前

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

    点赞 评论 复制链接分享
  • weixin_39698217 weixin_39698217 4月前

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

    https://github.com/pytroll/satpy/blob/ab17a717accfc0a7d789ab76912cd58ecad0dbd5/satpy/readers/grib.py#L145-L147

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

    点赞 评论 复制链接分享
  • weixin_39963255 weixin_39963255 4月前

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

    https://github.com/pytroll/satpy/blob/ab17a717accfc0a7d789ab76912cd58ecad0dbd5/satpy/readers/grib.py#L167

    and this part:

    https://github.com/pytroll/satpy/blob/ab17a717accfc0a7d789ab76912cd58ecad0dbd5/satpy/readers/grib.py#L185

    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:

    https://github.com/pytroll/satpy/blob/ab17a717accfc0a7d789ab76912cd58ecad0dbd5/satpy/readers/grib.py#L149-L150

    to:

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

    Your case might work.

    点赞 评论 复制链接分享
  • weixin_39698217 weixin_39698217 4月前

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

    点赞 评论 复制链接分享
  • weixin_39963255 weixin_39963255 4月前

    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.

    点赞 评论 复制链接分享
  • weixin_39698217 weixin_39698217 4月前

    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.

    点赞 评论 复制链接分享
  • weixin_39698217 weixin_39698217 4月前

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

    点赞 评论 复制链接分享
  • weixin_39698217 weixin_39698217 4月前

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

    python
    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?

    点赞 评论 复制链接分享
  • weixin_39963255 weixin_39963255 4月前

    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 https://github.com/pytroll/pyresample/issues/267 and point 1).

    点赞 评论 复制链接分享
  • weixin_39698217 weixin_39698217 4月前

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

    点赞 评论 复制链接分享
  • weixin_39963255 weixin_39963255 4月前

    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.

    点赞 评论 复制链接分享

相关推荐