pyresample icon indicating copy to clipboard operation
pyresample copied to clipboard

`get_sum` not matched with `bucket_sum`

Open zxdawn opened this issue 4 years ago • 7 comments

Code

import proplot as plot
import pyresample
import xarray as xr
from satpy import Scene
from pyresample import create_area_def
from pyresample.bucket import BucketResampler

abi_dir = '../data/GOES-16/ABI_L2/'
channel = 'HT'
reader = 'abi_l2_nc'

lon_min = -125
lon_max = -65
lat_min = 25
lat_max = 50
resample_res = 0.2 # resolution of resample grid (units: degree)

area = create_area_def('track_area',
                       {'proj': 'longlat', 'datum': 'WGS84'},
                       area_extent=[lon_min-resample_res/2,
                                    lat_min-resample_res/2,
                                    lon_max+resample_res/2,
                                    lat_max+resample_res/2],
                       resolution=resample_res,
                       units='degrees',
                       description=f'{resample_res}x{resample_res} degree lat-lon grid')

scn = Scene(['../data/GOES-16/ABI_L2/OR_ABI-L2-ACHAC-M6_G16_s20201530001131_e20201530003504_c20201530006018.nc'],
            reader=reader)
scn.load([channel])
new_scn = scn.resample(area, cache_dir=cache_dir, resample='bucket_sum')

lons, lats = scn[channel].area.get_lonlats()
resampler = BucketResampler(area,
                            xr.DataArray(lons).chunk(chunks={'dim_0':100}).data,
                            xr.DataArray(lats).chunk(chunks={'dim_0':100}).data)
sums = resampler.get_sum(scn[channel])

new_lons, new_lats = new_scn[channel].area.get_lonlats()
fig, axs = plot.subplots(axwidth=6, ncols=3, share=0)

scn[channel].plot(ax=axs[0], cmap='viridis')

m1 = axs[1].pcolormesh(new_lons, new_lats, new_scn[channel], cmap='viridis')
axs[1].colorbar(m1, loc='r')

m2 = axs[2].pcolormesh(new_lons, new_lats, sums, cmap='viridis')
axs[2].colorbar(m2, loc='r')

axs.format(collabels=['Original CTH', 'bucket_sum', 'Manual summation'])

fig.savefig('comparison.jpg')

Problem description

  1. Abnormal strips of the get_sum plot
  2. Missing data in the lower-left corner of the bucket_sum plot
  3. get_sum not matched with bucket_sum

Actual Result, Traceback if applicable

comparison

Versions of Python, package at hand and relevant dependencies

satpy 0.27.0
pyresample 1.19.0

zxdawn avatar May 08 '21 08:05 zxdawn

Took a while to figure this out.

new_scn = scn.resample(area, cache_dir=cache_dir, resample='bucket_sum')

has a typo. It should be resampler='bucket_sum'. With the typo, Satpy defaults to the 'nearest' resampler. The last image in the triplet is as it should.

pnuu avatar May 17 '21 08:05 pnuu

Oh, and the bucket resamplers don't have caching, so the cache_dir is unnecessary.

pnuu avatar May 17 '21 08:05 pnuu

@pnuu Thanks, but I got another error when using

new_scn = scn.resample(area, resampler='bucket_sum')
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-7-d8e0843844a5> in <module>
      2             reader=reader)
      3 scn.load([channel])
----> 4 new_scn = scn.resample(area,
      5 #                        cache_dir=cache_dir,
      6                        resampler='bucket_sum')

~/new/miniconda3/envs/pyresample_min/lib/python3.8/site-packages/satpy/scene.py in resample(self, destination, datasets, generate, unload, resampler, reduce_data, **resample_kwargs)
    819         # we may have some datasets we asked for but don't exist yet
    820         new_scn._wishlist = self._wishlist.copy()
--> 821         self._resampled_scene(new_scn, destination, resampler=resampler,
    822                               reduce_data=reduce_data, **resample_kwargs)
    823 

~/new/miniconda3/envs/pyresample_min/lib/python3.8/site-packages/satpy/scene.py in _resampled_scene(self, new_scn, destination_area, reduce_data, **resample_kwargs)
    775             kwargs = resample_kwargs.copy()
    776             kwargs['resampler'] = resamplers[source_area]
--> 777             res = resample_dataset(dataset, destination_area, **kwargs)
    778             new_datasets[ds_id] = res
    779             if ds_id in new_scn._datasets:

~/new/miniconda3/envs/pyresample_min/lib/python3.8/site-packages/satpy/resample.py in resample_dataset(dataset, destination_area, **kwargs)
   1376 
   1377     fill_value = kwargs.pop('fill_value', get_fill_value(dataset))
-> 1378     new_data = resample(source_area, dataset, destination_area, fill_value=fill_value, **kwargs)
   1379     new_attrs = new_data.attrs
   1380     new_data.attrs = dataset.attrs.copy()

~/new/miniconda3/envs/pyresample_min/lib/python3.8/site-packages/satpy/resample.py in resample(source_area, data, destination_area, resampler, **kwargs)
   1339         res = [resampler_instance.resample(ds, **kwargs) for ds in data]
   1340     else:
-> 1341         res = resampler_instance.resample(data, **kwargs)
   1342 
   1343     return res

~/new/miniconda3/envs/pyresample_min/lib/python3.8/site-packages/satpy/resample.py in resample(self, data, **kwargs)
   1118         else:
   1119             dims = data.dims
-> 1120         result = self.compute(data_arr, **kwargs)
   1121         coords = {}
   1122         if 'bands' in data.coords:

~/new/miniconda3/envs/pyresample_min/lib/python3.8/site-packages/satpy/resample.py in compute(self, data, skipna, **kwargs)
   1222                 results.append(res)
   1223         else:
-> 1224             res = self.resampler.get_sum(data, **kwargs)
   1225             results.append(res)
   1226 

TypeError: get_sum() got an unexpected keyword argument 'fill_value'

zxdawn avatar May 17 '21 08:05 zxdawn

Hi @zxdawn, thanks for spotting the bug, I must have introduced it with https://github.com/pytroll/satpy/pull/1539. 'fill_value' is not accepted atm by get_sum... before my PR, only specific kwargs were passed from satpy so fill_value was lost. Now all kwargs are passed and the bug is exposed. Not sure what is the best strategy to fix this, to avoid compatibility issues...

ameraner avatar May 17 '21 09:05 ameraner

Maybe satpy.resample._get_arg_to_pass_for_skipna_handling() could be adapted and renamed to filter out all unsupported kwargs?

pnuu avatar May 17 '21 09:05 pnuu

Another point: Shouldn't get_sum keep nan or fill_value instead of 0?

zxdawn avatar May 19 '21 08:05 zxdawn

My idea when doing the initial implementation was to have the real sum at each target location, thus I reasoned zero would be the correct value for the pixels where there is no data. This is arguable, but a fill_value kwarg defaulting to zero would be fine by me :-)

pnuu avatar May 28 '21 09:05 pnuu