openclimatefix / Satip

Satip contains the code necessary for retrieving, transforming and storing EUMETSAT data
https://satip.readthedocs.io/
MIT License
41 stars 29 forks source link

bug: infs in lats and lons #147

Open peterdudfield opened 1 year ago

peterdudfield commented 1 year ago

THere seem to be infs in the lats and lons

This occurred when the RSS data was back online on 20th December

<ipython-input-2-7bff5989e41b>:284: UserWarning: Attribute access to DataIDs is deprecated, use key access instead.
  logger.debug(f'Checking for infs in {scene._datasets.keys()[i].name}')
Traceback (most recent call last):
  File "<ipython-input-2-7bff5989e41b>", line 265, in convert_scene_to_dataarray
    scene = scene.crop(ll_bbox=GEOGRAPHIC_BOUNDS[area])
  File "/Users/peterdudfield/Documents/Github/nwp/venv/lib/python3.8/site-packages/satpy/scene.py", line 717, in crop
    new_coarsest_area, min_y_slice, min_x_slice = self._slice_area_from_bbox(
  File "/Users/peterdudfield/Documents/Github/nwp/venv/lib/python3.8/site-packages/satpy/scene.py", line 600, in _slice_area_from_bbox
    x_slice, y_slice = src_area.get_area_slices(dst_area)
  File "/Users/peterdudfield/Documents/Github/nwp/venv/lib/python3.8/site-packages/pyresample/geometry.py", line 2616, in get_area_slices
    intersection = data_boundary.contour_poly.intersection(
  File "/Users/peterdudfield/Documents/Github/nwp/venv/lib/python3.8/site-packages/pyresample/spherical.py", line 713, in intersection
    return self._bool_oper(other, -1)
  File "/Users/peterdudfield/Documents/Github/nwp/venv/lib/python3.8/site-packages/pyresample/spherical.py", line 649, in _bool_oper
    arcs1 = [edge for edge in self.aedges()]
  File "/Users/peterdudfield/Documents/Github/nwp/venv/lib/python3.8/site-packages/pyresample/spherical.py", line 649, in <listcomp>
    arcs1 = [edge for edge in self.aedges()]
  File "/Users/peterdudfield/Documents/Github/nwp/venv/lib/python3.8/site-packages/pyresample/spherical.py", line 571, in aedges
    yield Arc(SCoordinate(lon_start, lat_start),
  File "/Users/peterdudfield/Documents/Github/nwp/venv/lib/python3.8/site-packages/pyresample/spherical.py", line 131, in __init__
    lon, lat = _check_lon_lat(lon, lat)
  File "/Users/peterdudfield/Documents/Github/nwp/venv/lib/python3.8/site-packages/pyresample/spherical.py", line 119, in _check_lon_lat
    _check_lat_validity(lat)
  File "/Users/peterdudfield/Documents/Github/nwp/venv/lib/python3.8/site-packages/pyresample/spherical.py", line 109, in _check_lat_validity
    raise ValueError("Latitude values can not contain inf values.")
ValueError: Latitude values can not contain inf values.
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
  File "/Users/peterdudfield/Documents/Github/nwp/venv/lib/python3.8/site-packages/IPython/core/interactiveshell.py", line 3251, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-7-d4b6a84371db>", line 1, in <module>
    dataarray: xr.DataArray = convert_scene_to_dataarray(
  File "<ipython-input-2-7bff5989e41b>", line 293, in convert_scene_to_dataarray
    scene = scene.crop(ll_bbox=GEOGRAPHIC_BOUNDS[area])
  File "/Users/peterdudfield/Documents/Github/nwp/venv/lib/python3.8/site-packages/satpy/scene.py", line 717, in crop
    new_coarsest_area, min_y_slice, min_x_slice = self._slice_area_from_bbox(
  File "/Users/peterdudfield/Documents/Github/nwp/venv/lib/python3.8/site-packages/satpy/scene.py", line 600, in _slice_area_from_bbox
    x_slice, y_slice = src_area.get_area_slices(dst_area)
  File "/Users/peterdudfield/Documents/Github/nwp/venv/lib/python3.8/site-packages/pyresample/geometry.py", line 2616, in get_area_slices
    intersection = data_boundary.contour_poly.intersection(
  File "/Users/peterdudfield/Documents/Github/nwp/venv/lib/python3.8/site-packages/pyresample/spherical.py", line 713, in intersection
    return self._bool_oper(other, -1)
  File "/Users/peterdudfield/Documents/Github/nwp/venv/lib/python3.8/site-packages/pyresample/spherical.py", line 649, in _bool_oper
    arcs1 = [edge for edge in self.aedges()]
  File "/Users/peterdudfield/Documents/Github/nwp/venv/lib/python3.8/site-packages/pyresample/spherical.py", line 649, in <listcomp>
    arcs1 = [edge for edge in self.aedges()]
  File "/Users/peterdudfield/Documents/Github/nwp/venv/lib/python3.8/site-packages/pyresample/spherical.py", line 571, in aedges
    yield Arc(SCoordinate(lon_start, lat_start),
  File "/Users/peterdudfield/Documents/Github/nwp/venv/lib/python3.8/site-packages/pyresample/spherical.py", line 131, in __init__
    lon, lat = _check_lon_lat(lon, lat)
  File "/Users/peterdudfield/Documents/Github/nwp/venv/lib/python3.8/site-packages/pyresample/spherical.py", line 119, in _check_lon_lat
    _check_lat_validity(lat)
  File "/Users/peterdudfield/Documents/Github/nwp/venv/lib/python3.8/site-packages/pyresample/spherical.py", line 109, in _check_lat_validity
    raise ValueError("Latitude values can not contain inf values.")
ValueError: Latitude values can not contain inf values.
peterdudfield commented 1 year ago

I tried doing something like this, but it did not work

        except ValueError("Latitude values can not contain inf values."):
            """ This can sometimes happen where the area has a inf values with 'area_extent_ll'

            In general I (PD) dont like this, as 
            1. I think this should be sorted in satpy
            2. we are using _datasets which we shouldn't
            3. Seems dodgy just replacing 'inf' values.
            """
            n_datasets = len(scene._datasets)
            for i in range(n_datasets):
                logger.debug(f'Checking for infs in {scene._datasets.keys()[i].name}')
                ds = scene._datasets[scene._datasets.keys()[i]]
                area = ds.attrs.get('area')
                if area.area_extent_ll == (np.inf,np.inf,np.inf,np.inf):
                    logger.debug('Found inf values in "area_extent_ll", so changing them to "area_extent"')
                    area.area_extent_ll = area.area_extent

                    scene._datasets[scene._datasets.keys()[i]].attrs['area'] = area

            scene = scene.crop(ll_bbox=GEOGRAPHIC_BOUNDS[area])

this would be in satip/utils 273

peterdudfield commented 1 year ago

Idea: to have an easy way to switch to the 15 mins satellite data. I.e oput this in the click argument

peterdudfield commented 1 year ago

I downgraded to pyresample==1.25.0 and seemed to not get the error locally, will ahve to test of course

peterdudfield commented 1 year ago

This seemed to sovle it, Ill try and put a bug report on pyresample, but its good to find out what version jump did this. Reference @devsjc and @jacobbieker just incase your interested