Add utility functions for generating GeoViews plots

doc/source/multiscene.rst

@@ -13,12 +13,17 @@ examples will walk through some basic use cases of the MultiScene.
     These features are still early in development and may change overtime as
     more user feedback is received and more features added.
 
+Blending Scenes in MultiScene
+-----------------------------
+Scenes contained in a MultiScene can be combined in different ways.
+
 Stacking scenes
 ---------------
 
-The MultiScene makes it easy to take multiple Scenes and stack them on top of
-each other. The code below takes two separate orbits from a VIIRS sensor and
-stacks them on top of each other.
+The code below uses the `stack` function of the :meth:`~MultiScene.blend`
+method to stack two separate orbits from a VIIRS sensor. The `stack` function uses
+the first dataset as the base and then iteratively fills in missing values from
+all other datasets.
 
     >>> from satpy import Scene, MultiScene
     >>> from glob import glob
@@ -34,6 +39,13 @@ stacks them on top of each other.
     >>> blended_scene = new_mscn.blend()
     >>> blended_scene.save_datasets()
 
+Timeseries
+----------
+Using the :meth:`~MultiScene.blend` method with the `timeseries` function will combine              
+multiple scenes from different time slots by time. A single `Scene` with each dataset/channel       
+extended by the time dimension will be returned. If used together with the :meth:`~Scene.to_geoviews`
+method, creation of interactive timeseries Bokeh plots is possible.
+
 Saving frames of an animation
 -----------------------------
 

doc/source/quickstart.rst

@@ -2,8 +2,8 @@
 Quickstart
 ==========
 
-Loading data
-============
+Loading and accessing data
+==========================
 
 .. testsetup:: *
     >>> import sys
@@ -127,11 +127,31 @@ To access the loaded data use the wavelength or name:
 
     >>> print(global_scene[0.6])
 
-To visualize loaded data in a pop-up window:
-
-    >>> global_scene.show(0.6)
-
-To make combine datasets and make a new dataset:
+Visualizing data                                                                                    
+================                                                                                    
+
+To visualize loaded data in a pop-up window:                                                        
+
+    >>> global_scene.show(0.6)                                                                      
+
+Alternatively if working in a Jupyter notebook the scene can be converted to a `geoviews            
+<http://geo.holoviews.org/index.html>` object using the :meth:`~satpy.scene.Scene.to_geoviews`      
+method. The geoviews package is not a rquirement of the satpy install so in order to use this
+feature the user needs to install the geoviews package himself.
+
+    >>> import holoviews as hv                                                                      
+    >>> import geoviews as gv                                                                       
+    >>> import geoviews.feature as gf                                                               
+    >>> gv.extension("bokeh", "matplotlib")                                                         
+    >>> %opts QuadMesh Image [width=600 height=400 colorbar=True] Feature [apply_ranges=False]      
+    >>> %opts Image QuadMesh (cmap='RdBu_r')                                                        
+    >>> gview = global_scene.to_geoviews(vdims = [0.6])                                             
+    >>> gview[::5,::5] * gf.coastline * gf.borders                                                  
+
+Creating new datasets                                                                               
+=====================                                                                               
+
+Calculations based on loaded datasets/channels can easily be assigned to a new dataset:
 
     >>> global_scene["ndvi"] = (global_scene[0.8] - global_scene[0.6]) / (global_scene[0.8] + global_scene[0.6])
     >>> global_scene.show("ndvi")

satpy/multiscene.py

@@ -28,6 +28,7 @@
 import dask
 import dask.array as da
 import xarray as xr
+import pandas as pd
 from satpy.scene import Scene
 from satpy.writers import get_enhanced_image
 from itertools import chain
@@ -99,6 +100,17 @@ def stack(datasets):
     return base
 
 
+def timeseries(datasets):
+    """Expands dataset with and concats by time dimension"""
+    expanded_ds = []
+    for ds in datasets:
+        tmp = ds.expand_dims("time")
+        tmp.coords["time"] = pd.DatetimeIndex([ds.attrs["start_time"]])
+        expanded_ds.append(tmp)
+
+    return xr.concat(expanded_ds, dim="time")
+
+
 class _SceneGenerator(object):
     """Fancy way of caching Scenes from a generator."""
 

satpy/scene.py

@@ -30,13 +30,16 @@
 from satpy.composites import CompositorLoader, IncompatibleAreas
 from satpy.config import get_environ_config_dir
 from satpy.dataset import (DatasetID, MetadataObject, dataset_walker,
-                           replace_anc)
+                           replace_anc, combine_metadata)
 from satpy.node import DependencyTree
 from satpy.readers import DatasetDict, load_readers
 from satpy.resample import (resample_dataset,
                             prepare_resampler, get_area_def)
 from satpy.writers import load_writer
 from pyresample.geometry import AreaDefinition, BaseDefinition
+from pyresample.utils import proj4_str_to_dict
+
+import xarray as xr
 from xarray import DataArray
 import numpy as np
 import six
@@ -1024,6 +1027,60 @@ def show(self, dataset_id, overlay=None):
             img.show()
         return img
 
+    def to_geoviews(self, gvtype=None, datasets=None, kdims=None, vdims=None, dynamic=False):
+        """Convert satpy Scene to geoviews.
+
+        Parameters
+        ----------
+        gvtype : gv plot type, optional, default gv.Image
+            One of gv.Image, gv.LineContours, gv.FilledContours, gv.Points
+            See Geoviews documentation for details.
+        datasets : list of string
+        kdims : list of str, optional
+            Key dimensions. See geoviews documentation for more information.
+        vdims : list of str, optional
+            Value dimensions. See geoviews documentation for more information.
+            If not given defaults to first data variable
+        dynamic : boolean, optional, default False
+
+        Returns
+        -------
+        geoviews object
+
+        Todo
+        ----
+        - better handling of projection information in datasets which are
+          to be passed to geoviews
+        """
+        try:
+            import geoviews as gv
+            from cartopy import crs
+        except ImportError:
+            import warnings
+            warnings.warn("This method needs the geoviews package installed.")
+
+        if gvtype is None:
+            gvtype = gv.Image
+
+        ds = self.to_xarray_dataset(datasets)
+
+        if vdims is None:
+            # by default select first data variable as display variable
+            vdims = ds.data_vars[ds.data_vars.keys()[0]].name
+
+        if hasattr(ds, "area"):
+            dscrs = ds.area.to_cartopy_crs()
+            gvds = gv.Dataset(ds, crs=dscrs)
+        else:
+            LOG.error("No area found in dataset")
+
+        if "latitude" in ds.coords.keys():
+            gview = gvds.to(gv.QuadMesh, kdims=["longitude", "latitude"], vdims=vdims, dynamic=dynamic)
+        else:
+            gview = gvds.to(gvtype, kdims=["x", "y"], vdims=vdims, dynamic=dynamic)
+
+        return gview
+
     def images(self):
         """Generate images for all the datasets from the scene."""
         for ds_id, projectable in self.datasets.items():
@@ -1060,6 +1117,37 @@ def save_datasets(self, writer="geotiff", datasets=None, compute=True, **kwargs)
         writer, save_kwargs = load_writer(writer, ppp_config_dir=self.ppp_config_dir, **kwargs)
         return writer.save_datasets(datasets, compute=compute, **save_kwargs)
 
+    def to_xarray_dataset(self, datasets=None):
+        """Merge all xr.DataArrays of a scene to a xr.DataSet.
+
+        Parameters
+        ----------
+        datasets : list of string, optional
+            List of datasets to include in the xr.Dataset
+
+        Returns
+        -------
+        xr.DataSet
+        """
+        dslist = [(i, j) for i, j in self.datasets.items() if i.name not in ["latitude", "longitude"]]
+
+        if datasets is not None:
+            dslist = [(i, j) for i, j in dslist if i.name in datasets]
+
+        ds_dict = {i.name: it.rename(i.name) for (i, it) in dslist}
+        mdata = combine_metadata(*tuple(i.attrs for i in self.datasets.values()))
+
+        if "latitude" in [k.name for k in self.datasets.keys()]:
+            ds = xr.Dataset(ds_dict, coords={"latitude": (["y", "x"], self["latitude"]),
+                                             "longitude": (["y", "x"], self["longitude"],)})
+            # ds.longitude.values[ds.longitude.values>180] -= 360
+        else:
+            ds = xr.merge(ds_dict.values())
+
+        ds.attrs = mdata
+
+        return ds
+
     @classmethod
     def get_writer_by_ext(cls, extension):
         """Find the writer matching the *extension*."""

satpy/utils.py

@@ -288,3 +288,4 @@ def atmospheric_path_length_correction(data, cos_zen, limit=88., max_sza=95.):
     corr = corr.where(cos_zen.notnull(), 0)
 
     return data * corr
+