Resampling to a SwathDefinition fails

[jleinonen]

Describe the bug Resampling a scene to a pyresample.geometry.SwathDefinition fails if the lons and lats attributes to SwathDefinition have been passed as NumPy arrays. It works if they are passed as xarray.DataArray.

This behavior is inconsistent with the SwathDefinition docstring that specifies that lons and lats should be NumPy arrays.

Furthermore, once triggered to fail one time, it fails even with DataArray inputs until python is restarted (probably because of some caching I don't understand).

To Reproduce

# Your code here
from satpy import Scene
from pyresample.geometry import SwathDefinition
import numpy as np

files = [] # specify files here
scn = Scene(reader="abi_l1b", filenames=files) # I have GOES-R files so I use this reader
scn.load(scn.available_dataset_names())

lons = np.array([[-88.0, -87.0], [-88.0, -87.0]])
lats = np.array([[33.0, 33.0], [34.0, 34.0]])
swath = SwathDefinition(lons, lats)

local_scn = scn.resample(swath) # ** this throws an exception **

If we restart python first after triggering the above bug, the following works:

# Your code here
from satpy import Scene
from pyresample.geometry import SwathDefinition
import xarray

files = [] # specify files here
scn = Scene(reader="abi_l1b", filenames=files) # I have GOES-R files so I use this reader
scn.load(scn.available_dataset_names())

lons = xarray.DataArray([[-88.0, -87.0], [-88.0, -87.0]])
lats =  xarray.DataArray([[33.0, 33.0], [34.0, 34.0]])
swath = SwathDefinition(lons, lats)

local_scn = scn.resample(swath) # ** this works **

Expected behavior The scene scn should be resampled to the SwathDefinition (both examples above should work).

Actual results The first example above gives the following exception and traceback:

local_scn = scn.resample(swath)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-16-da2bd878a703> in <module>
----> 1 local_scn = scn.resample(swath)

<path>/lib/python3.8/site-packages/satpy/scene.py in resample(self, destination, datasets, generate, unload, resampler, reduce_data, **resample_kwargs)
   1120         # we may have some datasets we asked for but don't exist yet
   1121         new_scn._wishlist = self._wishlist.copy()
-> 1122         self._resampled_scene(new_scn, destination, resampler=resampler,
   1123                               reduce_data=reduce_data, **resample_kwargs)
   1124 

<path>/lib/python3.8/site-packages/satpy/scene.py in _resampled_scene(self, new_scn, destination_area, reduce_data, **resample_kwargs)
   1076             kwargs = resample_kwargs.copy()
   1077             kwargs['resampler'] = resamplers[source_area]
-> 1078             res = resample_dataset(dataset, destination_area, **kwargs)
   1079             new_datasets[ds_id] = res
   1080             if ds_id in new_scn._datasets:

<path>/lib/python3.8/site-packages/satpy/resample.py in resample_dataset(dataset, destination_area, **kwargs)
   1318 
   1319     fill_value = kwargs.pop('fill_value', get_fill_value(dataset))
-> 1320     new_data = resample(source_area, dataset, destination_area, fill_value=fill_value, **kwargs)
   1321     new_attrs = new_data.attrs
   1322     new_data.attrs = dataset.attrs.copy()

<path>/lib/python3.8/site-packages/satpy/resample.py in resample(source_area, data, destination_area, resampler, **kwargs)
   1281         res = [resampler_instance.resample(ds, **kwargs) for ds in data]
   1282     else:
-> 1283         res = resampler_instance.resample(data, **kwargs)
   1284 
   1285     return res

<path>/lib/python3.8/site-packages/satpy/resample.py in resample(self, data, cache_dir, mask_area, **kwargs)
    421 
    422         cache_id = self.precompute(cache_dir=cache_dir, **kwargs)
--> 423         return self.compute(data, cache_id=cache_id, **kwargs)
    424 
    425     def _create_cache_filename(self, cache_dir=None, prefix='',

<path>/lib/python3.8/site-packages/satpy/resample.py in compute(***failed resolving arguments***)
    604         LOG.debug("Resampling %s", str(data.name))
    605         res = self.resampler.get_sample_from_neighbour_info(data, fill_value)
--> 606         return update_resampled_coords(data, res, self.target_geo_def)
    607 
    608 

<path>/lib/python3.8/site-packages/satpy/resample.py in update_resampled_coords(old_data, new_data, new_area)
    333 
    334     # add crs, x, and y coordinates
--> 335     new_data = add_crs_xy_coords(new_data, new_area)
    336     return new_data
    337 

<path>/lib/python3.8/site-packages/satpy/resample.py in add_crs_xy_coords(data_arr, area)
    294         lons = area.lons
    295         lats = area.lats
--> 296         lons.attrs.setdefault('standard_name', 'longitude')
    297         lons.attrs.setdefault('long_name', 'longitude')
    298         lons.attrs.setdefault('units', 'degrees_east')

AttributeError: 'numpy.ndarray' object has no attribute 'attrs'

Workaround If you pass the lons and lats as DataArray, it seems to work fine.

Environment Info:

• OS: Linux (but probably not platform dependent)
• Satpy Version: 0.23.0
• PyResample Version: 1.16.0

[djhoese]

This is very interesting. This works with DataArrays because these are the type of SwathDefinitions that satpy produces. So I'm glad that that at least works. It looks like we need to do 2 things:

1. Update satpy to work with all SwathDefinitions or update the SwathDefinition class to provide these types of attributes when not backed by DataArrays.
2. Update the SwathDefinition docstring to talk about xarray and DataArrays not that it can support them (at least a little).

[wodowiesel]

i get the same problem as described too.

[gerritholl]

Resampling from a SwathDefinition also fails — but not exactly in the same place:

from satpy.resample import resample_dataset
from pyresample import create_area_def
from pyresample.geometry import SwathDefinition
import xarray as xr
import numpy as np
test_areadef = create_area_def(
        "test", 4087,
        area_extent=[0, 0, 500_000, 500_000],
        resolution=100_000)
test_swathdef = SwathDefinition(
        np.array([
            [0.4, 1.3, 2.2, 3.1, 4. ],
            [0.4, 1.3, 2.2, 3.1, 4. ],
            [0.4, 1.3, 2.2, 3.1, 4. ],
            [0.4, 1.3, 2.2, 3.1, 4. ],
            [0.4, 1.3, 2.2, 3.1, 4. ]]),
        np.array([
            [4. , 4. , 4. , 4. , 4. ],
            [3.1, 3.1, 3.1, 3.1, 3.1],
            [2.2, 2.2, 2.2, 2.2, 2.2],
            [1.3, 1.3, 1.3, 1.3, 1.3],
            [0.4, 0.4, 0.4, 0.4, 0.4]]))
ds = xr.DataArray(
        np.full((5, 5), 250),
        dims=("y", "x"),
        attrs={"area": test_swathdef})
new = resample_dataset(ds, test_areadef)

Result:

Traceback (most recent call last):
  File "/data/gholl/checkouts/protocode/resample-with-swathdef.py", line 27, in <module>
    new = resample_dataset(ds, test_areadef)
  File "/data/gholl/checkouts/satpy/satpy/resample.py", line 1395, in resample_dataset
    new_data = resample(source_area, dataset, destination_area, fill_value=fill_value, **kwargs)
  File "/data/gholl/checkouts/satpy/satpy/resample.py", line 1358, in resample
    res = resampler_instance.resample(data, **kwargs)
  File "/data/gholl/checkouts/satpy/satpy/resample.py", line 428, in resample
    geo_dims = self.source_geo_def.lons.dims
AttributeError: 'numpy.ndarray' object has no attribute 'dims'. Did you mean: 'dumps'?

[gerritholl]

When resampling in the other direction, just making lons/lats be xarray.DataArray is not enough. They must also have dimensions ('y', 'x') set explicitly. See #1997 for details.