pytroll / satpy

Python package for earth-observing satellite data processing
http://satpy.readthedocs.org/en/latest/
GNU General Public License v3.0
1.07k stars 295 forks source link

Return wrong latitude & longitude while reading FY-4B GHI data #2468

Open BigShuiTai opened 1 year ago

BigShuiTai commented 1 year ago

Describe the bug Return wrong latitude & longitude while reading FY-4B GHI data.

To Reproduce

# Your code here
files = [
    'FY4B-_GHI---_N_REGX_1330E_L1-_FDI-_MULT_NOM_20230415063800_20230415063859_2000M_V0001.HDF', 
    'FY4B-_GHI---_N_REGX_1330E_L1-_GEO-_MULT_NOM_20230415063800_20230415063859_2000M_V0001.HDF'
] # Files include L1 & GEO
scn = Scene(files, reader='ghi_l1')
scn.load(['C07'])
area_def = scn['C07'].attrs['area']
lons, lats = area_def.get_lonlats()
print(lons[(~np.isnan(lons))&(~np.isinf(lons))], lats[(~np.isnan(lats))&(~np.isinf(lats))])

Expected behavior Return correct lat/lon data that do not include any np.nan or np.inf.

Actual results

/usr/local/lib/python3.11/dist-packages/satpy/_config.py:129: DeprecationWarning: SelectableGroups dict interface is deprecated. Use select.
  for entry_point in cached_entry_points().get(name, []):
[DEBUG: 2023-05-06 22:45:56 : satpy.readers.yaml_reader] Reading ('/usr/local/lib/python3.11/dist-packages/satpy/etc/readers/ghi_l1.yaml',)
[DEBUG: 2023-05-06 22:45:56 : h5py._conv] Creating converter from 7 to 5
[DEBUG: 2023-05-06 22:45:56 : h5py._conv] Creating converter from 5 to 7
[DEBUG: 2023-05-06 22:45:56 : h5py._conv] Creating converter from 7 to 5
[DEBUG: 2023-05-06 22:45:56 : h5py._conv] Creating converter from 5 to 7
[DEBUG: 2023-05-06 22:45:56 : satpy.readers.yaml_reader] Assigning to ghi_l1: ['FY4B-_GHI---_N_REGX_1330E_L1-_FDI-_MULT_NOM_20230415063800_20230415063859_2000M_V0001.HDF', 'FY4B-_GHI---_N_REGX_1330E_L1-_GEO-_MULT_NOM_20230415063800_20230415063859_2000M_V0001.HDF']
[DEBUG: 2023-05-06 22:45:56 : satpy.composites.config_loader] Looking for composites config file ghi.yaml
[DEBUG: 2023-05-06 22:45:56 : satpy.composites.config_loader] Looking for composites config file visir.yaml
/usr/local/lib/python3.11/dist-packages/pkg_resources/__init__.py:121: DeprecationWarning: pkg_resources is deprecated as an API
  warnings.warn("pkg_resources is deprecated as an API", DeprecationWarning)
/usr/local/lib/python3.11/dist-packages/pkg_resources/__init__.py:2870: DeprecationWarning: Deprecated call to `pkg_resources.declare_namespace('mpl_toolkits')`.
Implementing implicit namespace packages (as specified in PEP 420) is preferred to `pkg_resources.declare_namespace`. See https://setuptools.pypa.io/en/latest/references/keywords.html#keyword-namespace-packages
  declare_namespace(pkg)
/usr/local/lib/python3.11/dist-packages/pkg_resources/__init__.py:2870: DeprecationWarning: Deprecated call to `pkg_resources.declare_namespace('mpl_toolkits.basemap_data')`.
Implementing implicit namespace packages (as specified in PEP 420) is preferred to `pkg_resources.declare_namespace`. See https://setuptools.pypa.io/en/latest/references/keywords.html#keyword-namespace-packages
  declare_namespace(pkg)
/usr/local/lib/python3.11/dist-packages/pkg_resources/__init__.py:2349: DeprecationWarning: Deprecated call to `pkg_resources.declare_namespace('mpl_toolkits')`.
Implementing implicit namespace packages (as specified in PEP 420) is preferred to `pkg_resources.declare_namespace`. See https://setuptools.pypa.io/en/latest/references/keywords.html#keyword-namespace-packages
  declare_namespace(parent)
/usr/local/lib/python3.11/dist-packages/pkg_resources/__init__.py:2870: DeprecationWarning: Deprecated call to `pkg_resources.declare_namespace('google')`.
Implementing implicit namespace packages (as specified in PEP 420) is preferred to `pkg_resources.declare_namespace`. See https://setuptools.pypa.io/en/latest/references/keywords.html#keyword-namespace-packages
  declare_namespace(pkg)
/usr/local/lib/python3.11/dist-packages/pkg_resources/__init__.py:2870: DeprecationWarning: Deprecated call to `pkg_resources.declare_namespace('zope')`.
Implementing implicit namespace packages (as specified in PEP 420) is preferred to `pkg_resources.declare_namespace`. See https://setuptools.pypa.io/en/latest/references/keywords.html#keyword-namespace-packages
  declare_namespace(pkg)
/usr/local/lib/python3.11/dist-packages/pkg_resources/__init__.py:2870: DeprecationWarning: Deprecated call to `pkg_resources.declare_namespace('ruamel')`.
Implementing implicit namespace packages (as specified in PEP 420) is preferred to `pkg_resources.declare_namespace`. See https://setuptools.pypa.io/en/latest/references/keywords.html#keyword-namespace-packages
  declare_namespace(pkg)
[DEBUG: 2023-05-06 22:45:56 : satpy.readers.ghi_l1] Reading in get_dataset C07.
[DEBUG: 2023-05-06 22:45:56 : satpy.readers.fy4_base] Calibrating to brightness_temperature
[DEBUG: 2023-05-06 22:45:56 : pyproj] PROJ_ERROR: geos: Invalid latitude
/usr/local/lib/python3.11/dist-packages/pyresample/geometry.py:1292: RuntimeWarning: invalid value encountered in multiply
  y = arange(*row_range, **y_kwargs) * -pixel_size_xy[1] + offset_xy[1]
[] []

Screenshots If applicable, add screenshots to help explain your problem.

Environment Info:

Additional context Add any other context about the problem here.

djhoese commented 1 year ago

I have no experience with this reader, but could you print the area_def you have and paste the output here?

BigShuiTai commented 1 year ago

I have no experience with this reader, but could you print the area_def you have and paste the output here?

Area ID: REGX_2000m
Description: FY-4 REGX area
Projection ID: FY-4, 2000m
Projection: {'a': '6378137.20703125', 'h': '35785862.7929687', 'lon_0': '133', 'no_defs': 'None', 'proj': 'geos', 'rf': '298.2561192803', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
Number of columns: 900
Number of rows: 1000
Area extent: (-3032011.4738, 2662009.9774, -1233004.5822, inf)
djhoese commented 1 year ago

Ah ok, there's our problem. The upper-most extent is inf. Not good.

The extent is calculated by converting:

        c_lats = self.file_content['/attr/Corner-Point Latitudes']
        c_lons = self.file_content['/attr/Corner-Point Longitudes']

To X/Y coordinates in the projection space and where Corner-Point Latitudes and longitudes are global attributes in the input files. Could you find these attributes in your files and tell me the values?

BigShuiTai commented 1 year ago

Ah ok, there's our problem. The upper-most extent is inf. Not good.

The extent is calculated by converting:

        c_lats = self.file_content['/attr/Corner-Point Latitudes']
        c_lons = self.file_content['/attr/Corner-Point Longitudes']

To X/Y coordinates in the projection space and where Corner-Point Latitudes and longitudes are global attributes in the input files. Could you find these attributes in your files and tell me the values?

The c_lats returns [6.55340000e+04 5.26822548e+01 2.61917515e+01 2.54671078e+01], and the c_lons returns [65534. 112.9593277 99.32941437 120.37347412]

I think 65534 is a wrong value for the GHI data.

djhoese commented 1 year ago

Oh yeah, that doesn't look good. @simonrp84 have you seen anything like this before?

simonrp84 commented 1 year ago

Yes, I have. The lats / lons produced by the reader are a fudge as the data from CMA didn't contain all the required information. It wouldn't surprise me if they've changed the contents of their data yet again and that's produced this problem.

Can take a look, but not for a couple of weeks.