pytroll / tutorial-satpy-half-day

Jupyter Notebook-based tutorial providing an overview of Satpy's features lasting about 4 hours
Other
55 stars 20 forks source link

Feedback on dask.persist, dynamic=True; Question about composite plot in matplotlib/cartopy #10

Open prasanjitdash opened 3 years ago

prasanjitdash commented 3 years ago

My set-up specs: Windows 10, 32 Gb RAM, 6 cores (12 logical), Python 3.8

Feedback/minor hiccups

  1. dask.persist throws tuple out of index error but commenting it out allowed to run the script successfully. https://github.com/pytroll/tutorial-satpy-half-day/blob/644a0b9419c628c00f66a94064d7cb70f11b10d6/notebooks/04_resampling.ipynb#L504

  2. dynamic=True created some issues (beyond my understanding!), but setting dynamic=False let it run. https://github.com/pytroll/tutorial-satpy-half-day/blob/644a0b9419c628c00f66a94064d7cb70f11b10d6/notebooks/07_cartopy_geoviews.ipynb#L246

Please take a look when possible. These are minor issues.

Question (critical) I have a question related to the plotting of a given composite using imshowrather than using xarray.plot. The [bands, row, col] order in the dask/xarray seems like an issue for imshowand needs to be in image order[row, col, bands].

For example, in the line: https://github.com/pytroll/tutorial-satpy-half-day/blob/644a0b9419c628c00f66a94064d7cb70f11b10d6/notebooks/05_composites.ipynb#L125

instead of using img_data.plot.imshow(rgb='bands', vmin=0, vmax=1)

I'd like to use cartopy+matplotlib as follows: cs = ax.imshow(img_data, transform=data_crs, vmin=0, vmax=1, zorder=2)

The img_data is a dask/xarray_core_array in [bands, row, col] order, e.g., (3 10000 10000). How to change it to image order [row, col, bands] e.g., (10000 10000 3)?

Overall, this is a terrific tool and very well implemented. Thank you, Prasanjit

prasanjitdash commented 3 years ago

https://github.com/pytroll/tutorial-satpy-half-day/blob/644a0b9419c628c00f66a94064d7cb70f11b10d6/notebooks/05_composites.ipynb#L125

instead of using img_data.plot.imshow(rgb='bands', vmin=0, vmax=1)

I'd like to use cartopy+matplotlib as follows: cs = ax.imshow(img_data, transform=data_crs, vmin=0, vmax=1, zorder=2)

Hi @djhoese

I think you can ignore the single question that I had. It seems like I can use CartoPy to plot a composite image. The single-band image was easy enough, but plotting the composite took me a while to understand. But if you know of a better way, I'd be glad to hear.

index

composite = 'true_color_lowres' new_scn is the resampled scene object. extent is the desired extent in [x0, x1, y0, y1]

%%time
with ProgressBar():
    target_crs = ccrs.PlateCarree()
    data_crs = new_scn[composite].attrs['area'].to_cartopy_crs()
    fig, ax = plt.subplots(nrows=1,ncols=1, subplot_kw={'projection': target_crs}, \
                       figsize=({insert fig size})
    ax.set_extent(extent)
    ax.add_feature(land, linewidth=1., linestyle='--', alpha=1, zorder=0);
    ax.add_feature(ocean, linewidth=1., linestyle='--', alpha=1, zorder=1);
    # cs = ax.imshow(img, transform=data_crs, vmin=0, vmax=1, zorder=2)
    img_data.plot.imshow(transform=data_crs, vmin=0, vmax=1, rgb='bands', zorder=2)
    gl = ax.gridlines(crs=target_crs, draw_labels=True, linewidth=0.5, color='gray', alpha=0.8, linestyle='--')
    gl.xformatter, gl.yformatter = LONGITUDE_FORMATTER, LATITUDE_FORMATTER
    ax.add_feature(coastline, linewidth=1., linestyle='--', alpha=1, zorder=3);
    plt.title(title, fontsize=18);
djhoese commented 3 years ago

As Johan mentioned on slack, there is this option:

from satpy.writers import get_enhanced_image
arr = np.moveaxis(get_enhanced_image(scn).data.values, 0, -1)
plt.imshow(arr, ...)

But I personally prefer letting xarray figure it out like you have it.

For your smaller issues, I'll try to check those out today. If you have the actual error messages produced that would help a lot. I know for the dynamic keyword argument there was a warning message we used to get about it being deprecated but we could never get it to work otherwise. I bet it is finally deprecated and functionality has changed. I'll see what I can figure out.

djhoese commented 3 years ago
  1. I am unable to reproduce any warning or error message with dask.persist. What's the error you're seeing?
  2. For geoviews, by default dynamic is set to False and I see this warning WARNING:param.Image04051: Setting non-parameter attribute dynamic=False using a mechanism intended only for parameters which is the warning I mentioned in my previous comment. Doing dynamic=True seemed to work fine for me. Could you describe in more detail the issues you're seeing?
prasanjitdash commented 3 years ago

As Johan mentioned on slack, there is this option:

from satpy.writers import get_enhanced_image
arr = np.moveaxis(get_enhanced_image(scn).data.values, 0, -1)
plt.imshow(arr, ...)

But I personally prefer letting xarray figure it out like you have it.

For your smaller issues, I'll try to check those out today. If you have the actual error messages produced that would help a lot. I know for the dynamic keyword argument there was a warning message we used to get about it being deprecated but we could never get it to work otherwise. I bet it is finally deprecated and functionality has changed. I'll see what I can figure out.

Thanks. as you rightly stated, xarray self managing it appears more appealing as it also takes the same arguments in. I am also happy to be aware of the alternative np.moveaxis option (just in case).

prasanjitdash commented 3 years ago
1. I am unable to reproduce any warning or error message with `dask.persist`. What's the error you're seeing?

2. For geoviews, by default dynamic is set to False and I see this warning `WARNING:param.Image04051: Setting non-parameter attribute dynamic=False using a mechanism intended only for parameters` which is the warning I mentioned in my previous comment. Doing `dynamic=True` seemed to work fine for me. Could you describe in more detail the issues you're seeing?

Thanks, again. These are minor issues and could be some version incompatibility as well but documenting it here in case someone else experiences the same.

  1. dask.persist error: (if I comment out those two lines, it works, of course defeats the intention to use the persist functionality)

IndexError Traceback (most recent call last)

in 4 with ProgressBar(): 5 abi_max_scn = abi_scn.resample(resampler='native') ----> 6 abi_max_scn['C05'], abi_max_scn['C07'] = dask.persist(abi_max_scn['C05'], abi_max_scn['C07']) 7 8 viirs_max_scn = viirs_scn.resample(abi_scn.max_area()) c:\users\dash\python38\lib\site-packages\dask\base.py in persist(*args, **kwargs) 774 results = schedule(dsk, keys, **kwargs) 775 d = dict(zip(keys, results)) --> 776 results2 = [r({k: d[k] for k in ks}, *s) for r, ks, s in postpersists] 777 return repack(results2) 778 c:\users\dash\python38\lib\site-packages\dask\base.py in (.0) 774 results = schedule(dsk, keys, **kwargs) 775 d = dict(zip(keys, results)) --> 776 results2 = [r({k: d[k] for k in ks}, *s) for r, ks, s in postpersists] 777 return repack(results2) 778 c:\users\dash\python38\lib\site-packages\xarray\core\dataarray.py in _dask_finalize(results, func, args, name) 777 @staticmethod 778 def _dask_finalize(results, func, args, name): --> 779 ds = func(results, *args) 780 variable = ds._variables.pop(_THIS_ARRAY) 781 coords = ds._variables c:\users\dash\python38\lib\site-packages\xarray\core\dataset.py in _dask_postpersist(dsk, info, *args) 781 if is_dask: 782 func, args2 = v --> 783 result = func(dsk, *args2) 784 else: 785 result = v c:\users\dash\python38\lib\site-packages\xarray\core\variable.py in _dask_finalize(results, array_func, array_args, dims, attrs, encoding) 436 def _dask_finalize(results, array_func, array_args, dims, attrs, encoding): 437 if isinstance(results, dict): # persist case --> 438 name = array_args[0] 439 results = {k: v for k, v in results.items() if k[0] == name} 440 data = array_func(results, *array_args) IndexError: tuple index out of range 2. `dynamic=True` error: (if I remove the keyword `dynamic`, then it works): --------------------------------------------------------------------------- KeyError Traceback (most recent call last) c:\users\dash\python38\lib\site-packages\IPython\core\formatters.py in __call__(self, obj, include, exclude) 968 969 if method is not None: --> 970 return method(include=include, exclude=exclude) 971 return None 972 else: c:\users\dash\python38\lib\site-packages\holoviews\core\dimension.py in _repr_mimebundle_(self, include, exclude) 1315 combined and returned. 1316 """ -> 1317 return Store.render(self) 1318 1319 c:\users\dash\python38\lib\site-packages\holoviews\core\options.py in render(cls, obj) 1403 data, metadata = {}, {} 1404 for hook in hooks: -> 1405 ret = hook(obj) 1406 if ret is None: 1407 continue c:\users\dash\python38\lib\site-packages\holoviews\ipython\display_hooks.py in pprint_display(obj) 280 if not ip.display_formatter.formatters['text/plain'].pprint: 281 return None --> 282 return display(obj, raw_output=True) 283 284 c:\users\dash\python38\lib\site-packages\holoviews\ipython\display_hooks.py in display(obj, raw_output, **kwargs) 256 elif isinstance(obj, (HoloMap, DynamicMap)): 257 with option_state(obj): --> 258 output = map_display(obj) 259 elif isinstance(obj, Plot): 260 output = render(obj) c:\users\dash\python38\lib\site-packages\holoviews\ipython\display_hooks.py in wrapped(element) 144 try: 145 max_frames = OutputSettings.options['max_frames'] --> 146 mimebundle = fn(element, max_frames=max_frames) 147 if mimebundle is None: 148 return {}, {} c:\users\dash\python38\lib\site-packages\holoviews\ipython\display_hooks.py in map_display(vmap, max_frames) 204 return None 205 --> 206 return render(vmap) 207 208 c:\users\dash\python38\lib\site-packages\holoviews\ipython\display_hooks.py in render(obj, **kwargs) 66 renderer = renderer.instance(fig='png') 67 ---> 68 return renderer.components(obj, **kwargs) 69 70 c:\users\dash\python38\lib\site-packages\holoviews\plotting\renderer.py in components(self, obj, fmt, comm, **kwargs) 408 doc = Document() 409 with config.set(embed=embed): --> 410 model = plot.layout._render_model(doc, comm) 411 if embed: 412 return render_model(model, comm) c:\users\dash\python38\lib\site-packages\panel\viewable.py in _render_model(self, doc, comm) 453 if comm is None: 454 comm = state._comm_manager.get_server_comm() --> 455 model = self.get_root(doc, comm) 456 457 if config.embed: c:\users\dash\python38\lib\site-packages\panel\viewable.py in get_root(self, doc, comm, preprocess) 510 """ 511 doc = init_doc(doc) --> 512 root = self._get_model(doc, comm=comm) 513 if preprocess: 514 self._preprocess(root) c:\users\dash\python38\lib\site-packages\panel\layout\base.py in _get_model(self, doc, root, parent, comm) 120 if root is None: 121 root = model --> 122 objects = self._get_objects(model, [], doc, root, comm) 123 props = dict(self._init_params(), objects=objects) 124 model.update(**self._process_param_change(props)) c:\users\dash\python38\lib\site-packages\panel\layout\base.py in _get_objects(self, model, old_objects, doc, root, comm) 110 else: 111 try: --> 112 child = pane._get_model(doc, root, model, comm) 113 except RerenderError: 114 return self._get_objects(model, current_objects[:i], doc, root, comm) c:\users\dash\python38\lib\site-packages\panel\pane\holoviews.py in _get_model(self, doc, root, parent, comm) 237 plot = self.object 238 else: --> 239 plot = self._render(doc, comm, root) 240 241 plot.pane = self c:\users\dash\python38\lib\site-packages\panel\pane\holoviews.py in _render(self, doc, comm, root) 304 kwargs['comm'] = comm 305 --> 306 return renderer.get_plot(self.object, **kwargs) 307 308 def _cleanup(self, root): c:\users\dash\python38\lib\site-packages\holoviews\plotting\bokeh\renderer.py in get_plot(self_or_cls, obj, doc, renderer, **kwargs) 71 combining the bokeh model with another plot. 72 """ ---> 73 plot = super(BokehRenderer, self_or_cls).get_plot(obj, doc, renderer, **kwargs) 74 if plot.document is None: 75 plot.document = Document() if self_or_cls.notebook_context else curdoc() c:\users\dash\python38\lib\site-packages\holoviews\plotting\renderer.py in get_plot(self_or_cls, obj, doc, renderer, comm, **kwargs) 218 219 # Initialize DynamicMaps with first data item --> 220 initialize_dynamic(obj) 221 222 if not renderer: c:\users\dash\python38\lib\site-packages\holoviews\plotting\util.py in initialize_dynamic(obj) 252 continue 253 if not len(dmap): --> 254 dmap[dmap._initial_key()] 255 256 c:\users\dash\python38\lib\site-packages\holoviews\core\spaces.py in __getitem__(self, key) 1339 # Not a cross product and nothing cached so compute element. 1340 if cache is not None: return cache -> 1341 val = self._execute_callback(*tuple_key) 1342 if data_slice: 1343 val = self._dataslice(val, data_slice) c:\users\dash\python38\lib\site-packages\holoviews\core\spaces.py in _execute_callback(self, *args) 1108 1109 with dynamicmap_memoization(self.callback, self.streams): -> 1110 retval = self.callback(*args, **kwargs) 1111 return self._style(retval) 1112 c:\users\dash\python38\lib\site-packages\holoviews\core\spaces.py in __call__(self, *args, **kwargs) 712 713 try: --> 714 ret = self.callable(*args, **kwargs) 715 except KeyError: 716 # KeyError is caught separately because it is used to signal c:\users\dash\python38\lib\site-packages\holoviews\core\spaces.py in dynamic_mul(*args, **kwargs) 291 if isinstance(self, DynamicMap): 292 def dynamic_mul(*args, **kwargs): --> 293 element = self[args] 294 if reverse: 295 return other * element c:\users\dash\python38\lib\site-packages\holoviews\core\spaces.py in __getitem__(self, key) 1339 # Not a cross product and nothing cached so compute element. 1340 if cache is not None: return cache -> 1341 val = self._execute_callback(*tuple_key) 1342 if data_slice: 1343 val = self._dataslice(val, data_slice) c:\users\dash\python38\lib\site-packages\holoviews\core\spaces.py in _execute_callback(self, *args) 1108 1109 with dynamicmap_memoization(self.callback, self.streams): -> 1110 retval = self.callback(*args, **kwargs) 1111 return self._style(retval) 1112 c:\users\dash\python38\lib\site-packages\holoviews\core\spaces.py in __call__(self, *args, **kwargs) 712 713 try: --> 714 ret = self.callable(*args, **kwargs) 715 except KeyError: 716 # KeyError is caught separately because it is used to signal c:\users\dash\python38\lib\site-packages\holoviews\core\spaces.py in dynamic_mul(*args, **kwargs) 291 if isinstance(self, DynamicMap): 292 def dynamic_mul(*args, **kwargs): --> 293 element = self[args] 294 if reverse: 295 return other * element c:\users\dash\python38\lib\site-packages\holoviews\core\spaces.py in __getitem__(self, key) 1339 # Not a cross product and nothing cached so compute element. 1340 if cache is not None: return cache -> 1341 val = self._execute_callback(*tuple_key) 1342 if data_slice: 1343 val = self._dataslice(val, data_slice) c:\users\dash\python38\lib\site-packages\holoviews\core\spaces.py in _execute_callback(self, *args) 1108 1109 with dynamicmap_memoization(self.callback, self.streams): -> 1110 retval = self.callback(*args, **kwargs) 1111 return self._style(retval) 1112 c:\users\dash\python38\lib\site-packages\holoviews\core\spaces.py in __call__(self, *args, **kwargs) 712 713 try: --> 714 ret = self.callable(*args, **kwargs) 715 except KeyError: 716 # KeyError is caught separately because it is used to signal c:\users\dash\python38\lib\site-packages\holoviews\util\__init__.py in dynamic_operation(*key, **kwargs) 1040 1041 def dynamic_operation(*key, **kwargs): -> 1042 key, obj = resolve(key, kwargs) 1043 return apply(obj, *key, **kwargs) 1044 c:\users\dash\python38\lib\site-packages\holoviews\util\__init__.py in resolve(key, kwargs) 1029 elif isinstance(map_obj, DynamicMap) and map_obj._posarg_keys and not key: 1030 key = tuple(kwargs[k] for k in map_obj._posarg_keys) -> 1031 return key, map_obj[key] 1032 1033 def apply(element, *key, **kwargs): c:\users\dash\python38\lib\site-packages\holoviews\core\spaces.py in __getitem__(self, key) 1339 # Not a cross product and nothing cached so compute element. 1340 if cache is not None: return cache -> 1341 val = self._execute_callback(*tuple_key) 1342 if data_slice: 1343 val = self._dataslice(val, data_slice) c:\users\dash\python38\lib\site-packages\holoviews\core\spaces.py in _execute_callback(self, *args) 1108 1109 with dynamicmap_memoization(self.callback, self.streams): -> 1110 retval = self.callback(*args, **kwargs) 1111 return self._style(retval) 1112 c:\users\dash\python38\lib\site-packages\holoviews\core\spaces.py in __call__(self, *args, **kwargs) 712 713 try: --> 714 ret = self.callable(*args, **kwargs) 715 except KeyError: 716 # KeyError is caught separately because it is used to signal c:\users\dash\python38\lib\site-packages\holoviews\core\data\__init__.py in load_subset(*args) 990 def load_subset(*args): 991 constraint = dict(zip(dim_names, args)) --> 992 group = self.select(**constraint) 993 if np.isscalar(group): 994 return group_type(([group],), group=self.group, c:\users\dash\python38\lib\site-packages\holoviews\core\data\__init__.py in pipelined_fn(*args, **kwargs) 203 204 try: --> 205 result = method_fn(*args, **kwargs) 206 if PipelineMeta.disable: 207 return result c:\users\dash\python38\lib\site-packages\holoviews\core\data\__init__.py in select(self, selection_expr, selection_specs, **selection) 631 # Handle selection kwargs 632 if selection: --> 633 data = self.interface.select(self, **selection) 634 else: 635 data = self.data c:\users\dash\python38\lib\site-packages\holoviews\core\data\xarray.py in select(cls, dataset, selection_mask, **selection) 565 else: 566 validated[dim] = v --> 567 data = dataset.data.sel(**validated) 568 569 # Restore constant dimensions c:\users\dash\python38\lib\site-packages\xarray\core\dataset.py in sel(self, indexers, method, tolerance, drop, **indexers_kwargs) 2063 """ 2064 indexers = either_dict_or_kwargs(indexers, indexers_kwargs, "sel") -> 2065 pos_indexers, new_indexes = remap_label_indexers( 2066 self, indexers=indexers, method=method, tolerance=tolerance 2067 ) c:\users\dash\python38\lib\site-packages\xarray\core\coordinates.py in remap_label_indexers(obj, indexers, method, tolerance, **indexers_kwargs) 394 } 395 --> 396 pos_indexers, new_indexes = indexing.remap_label_indexers( 397 obj, v_indexers, method=method, tolerance=tolerance 398 ) c:\users\dash\python38\lib\site-packages\xarray\core\indexing.py in remap_label_indexers(data_obj, indexers, method, tolerance) 268 coords_dtype = data_obj.coords[dim].dtype 269 label = maybe_cast_to_coords_dtype(label, coords_dtype) --> 270 idxr, new_idx = convert_label_indexer(index, label, dim, method, tolerance) 271 pos_indexers[dim] = idxr 272 if new_idx is not None: c:\users\dash\python38\lib\site-packages\xarray\core\indexing.py in convert_label_indexer(index, label, index_name, method, tolerance) 187 indexer = index.get_loc(label.item()) 188 else: --> 189 indexer = index.get_loc( 190 label.item(), method=method, tolerance=tolerance 191 ) c:\users\dash\python38\lib\site-packages\pandas\core\indexes\datetimes.py in get_loc(self, key, method, tolerance) 681 else: 682 # unrecognized type --> 683 raise KeyError(key) 684 685 try: KeyError: 1526072428800000000
djhoese commented 3 years ago
  1. Regarding the resampling/dask error, what version of pyresample are you running?
  2. For geoviews, I'm really not sure. I wasn't able to reproduce this and that error seems very very strange (just a number?). Maybe we just say this is a version incompatibility thing. Satpy isn't doing anything too fancy to use these so maybe it is a pandas + xarray + geoviews/holoviews incompatibility that will get fixed if you upgrade or in the future when someone runs into something similar. There is a little bit too much involved with Satpy and everything to be able to file a bug with them unless we can narrow this down to just an xarray + geoviews simple example.
prasanjitdash commented 3 years ago
1. Regarding the resampling/dask error, what version of pyresample are you running?

2. For geoviews, I'm really not sure. I wasn't able to reproduce this and that error seems very very strange (just a number?). Maybe we just say this is a version incompatibility thing. Satpy isn't doing anything too fancy to use these so maybe it is a pandas + xarray + geoviews/holoviews incompatibility that will get fixed if you upgrade or in the future when someone runs into something similar. There is a little bit too much involved with Satpy and everything to be able to file a bug with them unless we can narrow this down to just an xarray + geoviews simple example.

Thanks, below are the versions of the Python stack. Totally agree; with just a number, it is hard to find out the root cause, and perhaps it is due to some incompatibility. Since this is low-hanging fruit, let's revisit this in some time, perhaps when we hear of similar issues. This reporting here is just for the record.

pyresample        1.20.0
satpy             0.29.0
pandas            1.2.3
xarray            0.15.1
geoviews          1.9.1
holoviews         1.14.5
dask              2021.7.1
python            3.8.0