Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

get_sum not matched with bucket_sum #354

Open
zxdawn opened this issue May 8, 2021 · 7 comments · May be fixed by #364
Open

get_sum not matched with bucket_sum #354

zxdawn opened this issue May 8, 2021 · 7 comments · May be fixed by #364
Assignees
@ameraner
Labels
bug

Comments

@zxdawn
Copy link
Member

zxdawn commented May 8, 2021 •
edited

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
@pnuu
Copy link
Member

pnuu commented May 17, 2021

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
Copy link
Member

pnuu commented May 17, 2021

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

@zxdawn
Copy link
Member Author

zxdawn commented May 17, 2021

@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'

@ameraner
Copy link
Member

Hi @zxdawn,
thanks for spotting the bug, I must have introduced it with pytroll/satpy#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...

@pnuu
Copy link
Member

pnuu commented May 17, 2021

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

@pnuu pnuu added the bug label May 17, 2021
@zxdawn
Copy link
Member Author

zxdawn commented May 19, 2021

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

@pnuu
Copy link
Member

pnuu commented May 28, 2021

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

@zxdawn zxdawn linked a pull request May 28, 2021 that will close this issue
Draft
5 tasks
@zxdawn zxdawn mentioned this issue Jun 30, 2021
Open
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees

@ameraner ameraner

Labels
bug
Projects
None yet
Milestone
No milestone
Development

Successfully merging a pull request may close this issue.

Set summed bin without data to fill_value zxdawn/pyresample
3 participants
@pnuu @zxdawn @ameraner