Add utility functions for generating GeoViews plots

.travis.yml

@@ -5,7 +5,7 @@ env:
         - PYTHON_VERSION=$TRAVIS_PYTHON_VERSION
         - NUMPY_VERSION=stable
         - MAIN_CMD='python setup.py'
-        - CONDA_DEPENDENCIES='xarray dask toolz Cython sphinx cartopy pillow matplotlib scipy pyyaml pyproj pyresample coveralls coverage codecov behave netcdf4 h5py h5netcdf gdal rasterio<1.0.14 imageio pyhdf mock libtiff pycoast pydecorate'
+        - CONDA_DEPENDENCIES='xarray dask toolz Cython sphinx cartopy pillow matplotlib scipy pyyaml pyproj pyresample coveralls coverage codecov behave netcdf4 h5py h5netcdf gdal rasterio<1.0.14 imageio pyhdf mock libtiff pycoast pydecorate geoviews'
         - PIP_DEPENDENCIES='trollsift trollimage pyspectral pyorbital libtiff'
         - SETUP_XVFB=False
         - EVENT_TYPE='push pull_request'

doc/source/conf.py

@@ -248,4 +248,5 @@ def __getattr__(cls, name):
     'trollsift': ('https://trollsift.readthedocs.io/en/stable', None),
     'trollimage': ('https://trollimage.readthedocs.io/en/stable', None),
     'pydecorate': ('https://pydecorate.readthedocs.io/en/stable', None),
+    'geoviews': ('http://geo.holoviews.org', None),
 }

doc/source/multiscene.rst

@@ -13,12 +13,18 @@ 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 :meth:`~satpy.multiscene.MultiScene.blend` method of
+the ``MultiScene`` object to stack two separate orbits from a VIIRS sensor. By
+default the ``blend`` method will use the :func:`~satpy.multiscene.stack`
+function which uses the first dataset as the base of the image and then
+iteratively overlays the remaining datasets on top.
 
     >>> from satpy import Scene, MultiScene
     >>> from glob import glob
@@ -34,6 +40,36 @@ stacks them on top of each other.
     >>> blended_scene = new_mscn.blend()
     >>> blended_scene.save_datasets()
 
+Timeseries
+**********
+
+Using the :meth:`~satpy.multiscene.MultiScene.blend` method with the
+:func:`~satpy.multiscene.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:`~satpy.scene.Scene.to_geoviews` method, creation of
+interactive timeseries Bokeh plots is possible.
+
+    >>> from satpy import Scene, MultiScene
+    >>> from satpy.multiscene import timeseries
+    >>> from glob import glob
+    >>> from pyresample.geometry import AreaDefinition
+    >>> my_area = AreaDefinition(...)
+    >>> scenes = [
+    ...    Scene(reader='viirs_sdr', filenames=glob('/data/viirs/day_1/*t180*.h5')),
+    ...    Scene(reader='viirs_sdr', filenames=glob('/data/viirs/day_2/*t180*.h5'))
+    ... ]
+    >>> mscn = MultiScene(scenes)
+    >>> mscn.load(['I04'])
+    >>> new_mscn = mscn.resample(my_area)
+    >>> blended_scene = new_mscn.blend(blend_function=timeseries)
+    >>> blended_scene['I04']
+    <xarray.DataArray (time: 2, y: 1536, x: 6400)>
+    dask.array<shape=(2, 1536, 6400), dtype=float64, chunksize=(1, 1536, 4096)>
+    Coordinates:
+      * time     (time) datetime64[ns] 2012-02-25T18:01:24.570942 2012-02-25T18:02:49.975797
+    Dimensions without coordinates: y, x
+
 Saving frames of an animation
 -----------------------------
 

doc/source/quickstart.rst

@@ -2,8 +2,8 @@
 Quickstart
 ==========
 
-Loading data
-============
+Loading and accessing data
+==========================
 
 .. testsetup:: *
     >>> import sys
@@ -106,7 +106,7 @@ SatPy allows loading file data by wavelengths in micrometers (shown above) or by
     >>> global_scene.load(["VIS006", "VIS008", "IR_108"])
 
 To have a look at the available channels for loading from your :class:`~satpy.scene.Scene` object use the
-:meth:`available_datasets <satpy.scene.Scene.available_dataset_names>` method:
+:meth:`~satpy.scene.Scene.available_dataset_names` method:
 
     >>> global_scene.available_dataset_names()
     ['HRV',
@@ -127,11 +127,32 @@ 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
+requirement of the base 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,8 +28,10 @@
 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 satpy.dataset import combine_metadata
 from itertools import chain
 
 try:
@@ -106,6 +108,19 @@ 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)
+
+    res = xr.concat(expanded_ds, dim="time")
+    res.attrs = combine_metadata(*[x.attrs for x in expanded_ds])
+    return res
+
+
 class _SceneGenerator(object):
     """Fancy way of caching Scenes from a generator."""
 

satpy/scene.py

@@ -30,13 +30,15 @@
 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.geometry import AreaDefinition, BaseDefinition, SwathDefinition
+
+import xarray as xr
 from xarray import DataArray
 import numpy as np
 import six
@@ -524,7 +526,13 @@ def slice(self, key):
         new_scn = self.copy()
         new_scn.wishlist = self.wishlist
         for area, dataset_ids in self.iter_by_area():
-            new_area = area[key] if area is not None else None
+            if area is not None:
+                # assume dimensions for area are y and x
+                one_ds = self[dataset_ids[0]]
+                area_key = tuple(sl for dim, sl in zip(one_ds.dims, key) if dim in ['y', 'x'])
+                new_area = area[area_key]
+            else:
+                new_area = None
             new_scn._slice_datasets(dataset_ids, key, new_area)
         return new_scn
 
@@ -1026,6 +1034,93 @@ 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.
+
+        Args:
+            gvtype (gv plot type):
+                One of gv.Image, gv.LineContours, gv.FilledContours, gv.Points
+                Default to :class:`geoviews.Image`.
+                See Geoviews documentation for details.
+            datasets (list): Limit included products to these datasets
+            kdims (list of str):
+                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  # noqa
+        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[list(ds.data_vars.keys())[0]].name
+
+        if hasattr(ds, "area") and hasattr(ds.area, 'to_cartopy_crs'):
+            dscrs = ds.area.to_cartopy_crs()
+            gvds = gv.Dataset(ds, crs=dscrs)
+        else:
+            gvds = gv.Dataset(ds)
+
+        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 to_xarray_dataset(self, datasets=None):
+        """Merge all xr.DataArrays of a scene to a xr.DataSet.
+
+        Parameters:
+            datasets (list):
+                List of products to include in the :class:`xarray.Dataset`
+
+        Returns: :class:`xarray.Dataset`
+
+        """
+        if datasets is not None:
+            datasets = [self[ds] for ds in datasets]
+        else:
+            datasets = [self.datasets.get(ds) for ds in self.wishlist]
+            datasets = [ds for ds in datasets if ds is not None]
+
+        ds_dict = {i.attrs['name']: i.rename(i.attrs['name']) for i in datasets if i.attrs.get('area') is not None}
+        mdata = combine_metadata(*tuple(i.attrs for i in datasets))
+        if mdata.get('area') is None or not isinstance(mdata['area'], SwathDefinition):
+            # either don't know what the area is or we have an AreaDefinition
+            ds = xr.merge(ds_dict.values())
+        else:
+            # we have a swath definition and should use lon/lat values
+            lons, lats = mdata['area'].get_lonlats()
+            if not isinstance(lons, DataArray):
+                lons = DataArray(lons, dims=('y', 'x'))
+                lats = DataArray(lats, dims=('y', 'x'))
+            # ds_dict['longitude'] = lons
+            # ds_dict['latitude'] = lats
+            ds = xr.Dataset(ds_dict, coords={"latitude": (["y", "x"], lats),
+                                             "longitude": (["y", "x"], lons)})
+
+        ds.attrs = mdata
+        return ds
+
     def images(self):
         """Generate images for all the datasets from the scene."""
         for ds_id, projectable in self.datasets.items():

satpy/tests/test_multiscene.py

@@ -190,12 +190,47 @@ def test_save_mp4(self):
                          'test_save_mp4_ds3_20180102_00_20180102_12.mp4')
 
 
+class TestBlendFuncs(unittest.TestCase):
+    """Test individual functions used for blending."""
+
+    def setUp(self):
+        """Set up test data."""
+        import xarray as xr
+        import dask.array as da
+        from datetime import datetime
+        from pyresample.geometry import AreaDefinition
+        area = AreaDefinition('test', 'test', 'test',
+                              {'proj': 'geos', 'lon_0': -95.5, 'h': 35786023.0},
+                              2, 2, [-200, -200, 200, 200])
+        ds1 = xr.DataArray(da.zeros((2, 2), chunks=-1), dims=('y', 'x'),
+                           attrs={'start_time': datetime(2018, 1, 1, 0, 0, 0), 'area': area})
+        self.ds1 = ds1
+        ds2 = xr.DataArray(da.zeros((2, 2), chunks=-1), dims=('y', 'x'),
+                           attrs={'start_time': datetime(2018, 1, 1, 1, 0, 0), 'area': area})
+        self.ds2 = ds2
+
+    def test_stack(self):
+        """Test the 'stack' function."""
+        from satpy.multiscene import stack
+        res = stack([self.ds1, self.ds2])
+        self.assertTupleEqual(self.ds1.shape, res.shape)
+
+    def test_timeseries(self):
+        """Test the 'timeseries' function."""
+        from satpy.multiscene import timeseries
+        import xarray as xr
+        res = timeseries([self.ds1, self.ds2])
+        self.assertIsInstance(res, xr.DataArray)
+        self.assertTupleEqual((2, self.ds1.shape[0], self.ds1.shape[1]), res.shape)
+
+
 def suite():
     """The test suite for test_multiscene."""
     loader = unittest.TestLoader()
     mysuite = unittest.TestSuite()
     mysuite.addTest(loader.loadTestsFromTestCase(TestMultiScene))
     mysuite.addTest(loader.loadTestsFromTestCase(TestMultiSceneSave))
+    mysuite.addTest(loader.loadTestsFromTestCase(TestBlendFuncs))
 
     return mysuite
 

satpy/tests/test_scene.py

@@ -1696,6 +1696,46 @@ def test_save_dataset_default(self):
             os.path.join(self.base_dir, 'test_20180101_000000.tif')))
 
 
+class TestSceneConversions(unittest.TestCase):
+    """Test Scene conversion to geoviews, xarray, etc."""
+
+    def test_geoviews_basic_with_area(self):
+        """Test converting a Scene to geoviews with an AreaDefinition."""
+        from satpy import Scene
+        import xarray as xr
+        import dask.array as da
+        from datetime import datetime
+        from pyresample.geometry import AreaDefinition
+        scn = Scene()
+        area = AreaDefinition('test', 'test', 'test',
+                              {'proj': 'geos', 'lon_0': -95.5, 'h': 35786023.0},
+                              2, 2, [-200, -200, 200, 200])
+        scn['ds1'] = xr.DataArray(da.zeros((2, 2), chunks=-1), dims=('y', 'x'),
+                                  attrs={'start_time': datetime(2018, 1, 1),
+                                         'area': area})
+        gv_obj = scn.to_geoviews()
+        # we assume that if we got something back, geoviews can use it
+        self.assertIsNotNone(gv_obj)
+
+    def test_geoviews_basic_with_swath(self):
+        """Test converting a Scene to geoviews with a SwathDefinition."""
+        from satpy import Scene
+        import xarray as xr
+        import dask.array as da
+        from datetime import datetime
+        from pyresample.geometry import SwathDefinition
+        scn = Scene()
+        lons = xr.DataArray(da.zeros((2, 2)))
+        lats = xr.DataArray(da.zeros((2, 2)))
+        area = SwathDefinition(lons, lats)
+        scn['ds1'] = xr.DataArray(da.zeros((2, 2), chunks=-1), dims=('y', 'x'),
+                                  attrs={'start_time': datetime(2018, 1, 1),
+                                         'area': area})
+        gv_obj = scn.to_geoviews()
+        # we assume that if we got something back, geoviews can use it
+        self.assertIsNotNone(gv_obj)
+
+
 def suite():
     """The test suite for test_scene."""
     loader = unittest.TestLoader()
@@ -1704,6 +1744,7 @@ def suite():
     mysuite.addTest(loader.loadTestsFromTestCase(TestSceneLoading))
     mysuite.addTest(loader.loadTestsFromTestCase(TestSceneResampling))
     mysuite.addTest(loader.loadTestsFromTestCase(TestSceneSaving))
+    mysuite.addTest(loader.loadTestsFromTestCase(TestSceneConversions))
 
     return mysuite