Introduction to DEA Coastlines notebook (#698)

* Add DEA Coastlines notebook and code, as well as metadata func

* Fix invalid html tags issue

* Run notebook from top

* Add notebook to readme

* Update existing coastal notebooks to better reflect DEA Coastlines code

* Minor markdown update

* Remove metadata load script for time being

* Re-run all cells

* Fix duplicated study area in intertidal notebook

* Update for review comments

* Remove try/except and add note about scikit-image issue

* Add additonal markdown and instructions for exporting data

* Add screenshots

* Reduce size of screenshots

DEA_datasets/README.rst

@@ -8,6 +8,7 @@ Notebooks introducing DEA's satellite datasets and derived products, including h
    :caption: DEA Datasets
 
    DEA_High_and_Low_Tide_Imagery.ipynb
+   DEA_Coastlines.ipynb
+   DEA_Waterbodies.ipynb
    Sentinel_2.ipynb
    TSmask.ipynb
-   DEA_Waterbodies.ipynb

Scripts/dea_plotting.py

@@ -29,12 +29,14 @@
 
 # Import required packages
 import math
+import branca
 import folium
 import calendar
 import ipywidgets
 import numpy as np
 import geopandas as gpd
 import matplotlib as mpl
+import matplotlib.cm as cm
 import matplotlib.patheffects as PathEffects
 import matplotlib.pyplot as plt
 import matplotlib.animation as animation
@@ -46,7 +48,6 @@
 from mpl_toolkits.axes_grid1 import make_axes_locatable
 from ipyleaflet import Map, Marker, Popup, GeoJSON, basemaps, Choropleth
 from skimage import exposure
-from branca.colormap import linear
 from odc.ui import image_aspect
 
 from matplotlib.animation import FuncAnimation
@@ -357,7 +358,7 @@ def display_map(x, y, crs='EPSG:4326', margin=-0.5, zoom_bias=0):
 def map_shapefile(gdf,
                   attribute,
                   continuous=False,
-                  colormap='YlOrRd_09',
+                  cmap='viridis',
                   basemap=basemaps.Esri.WorldImagery,
                   default_zoom=None,
                   hover_col=True,
@@ -382,13 +383,11 @@ def map_shapefile(gdf,
         Whether to plot data as a categorical or continuous variable. 
         Defaults to remapping the attribute which is suitable for 
         categorical data. For continuous data set `continuous` to True.
-    colormap : string, optional
-        Either a string giving the name of a `branca.colormap.linear` 
-        colormap or a `branca.colormap` object (for example, 
-        `branca.colormap.linear.YlOrRd_09`) that will be used to style 
-        the features in the GeoDataFrame. Features will be coloured 
-        based on the selected attribute. Defaults to the 'YlOrRd_09' 
-        colormap.
+    cmap : string, optional
+        A string giving the name of a `matplotlib.cm` colormap that will
+        be used to style the features in the GeoDataFrame. Features will
+        be coloured based on the selected attribute. Defaults to the 
+        'viridis' colormap.
     basemap : ipyleaflet.basemaps object, optional
         An optional `ipyleaflet.basemaps` object used as the basemap for 
         the interactive plot. Defaults to `basemaps.Esri.WorldImagery`.
@@ -403,7 +402,7 @@ def map_shapefile(gdf,
         custom shapefile field can be specified by supplying a string
         giving the name of the field to print. Set to False to prevent 
         any attributes from being printed.
-    **choropleth_kwargs :
+    **style_kwargs :
         Optional keyword arguments to pass to the `style` paramemter of
         the `ipyleaflet.Choropleth` function. This can be used to 
         control the appearance of the shapefile, for example 'stroke' 
@@ -484,15 +483,16 @@ def on_hover(event, id, properties):
     # overwritten/customised by `choropleth_kwargs` values
     style_kwargs = dict({'fillOpacity': 0.8}, **style_kwargs)
 
-    # Get colormap from either string or `branca.colormap` object
-    if type(colormap) == str:
-        colormap = getattr(linear, colormap)
+    # Get `branca.colormap` object from matplotlib string
+    cm_cmap = cm.get_cmap(cmap, 30)
+    colormap = branca.colormap.LinearColormap([cm_cmap(i) for 
+                                               i in np.linspace(0, 1, 30)])
 
     # Create the choropleth
     choropleth = Choropleth(geo_data=data_geojson,
                             choro_data=id_class_dict,
                             colormap=colormap,
-                            style={**style_kwargs})
+                            style=style_kwargs)
 
     # If the vector data contains line features, they will not be 
     # be coloured by default. To resolve this, we need to manually copy
@@ -509,7 +509,8 @@ def on_hover(event, id, properties):
             choropleth.data['features'][i]['properties']['style']['fillColor']
 
         # Add GeoJSON layer to map
-        feature_layer = GeoJSON(data=choropleth.data)
+        feature_layer = GeoJSON(data=choropleth.data,
+                                style=style_kwargs)
         m.add_layer(feature_layer)
 
     else:

Scripts/dea_spatialtools.py

@@ -26,7 +26,7 @@
     largest_region
     transform_geojson_wgs_to_epsg
 
-Last modified: June 2020
+Last modified: November 2020
 
 '''
 
@@ -340,7 +340,7 @@ def subpixel_contours(da,
     `attribute_df` parameter can be used to pass custom attributes 
     to the output contour features.
     
-    Last modified: June 2020
+    Last modified: November 2020
     
     Parameters
     ----------  
@@ -418,10 +418,12 @@ def contours_to_multiline(da_i, z_value, min_vertices=2):
         '''
 
         # Extracts contours from array, and converts each discrete
-        # contour into a Shapely LineString feature
+        # contour into a Shapely LineString feature. If the function 
+        # returns a KeyError, this may be due to an unresolved issue in
+        # scikit-image: https://github.com/scikit-image/scikit-image/issues/4830
         line_features = [LineString(i[:,[1, 0]]) 
-                         for i in find_contours(da_i.data, z_value)
-                         if i.shape[0] > min_vertices]
+                         for i in find_contours(da_i.data, z_value) 
+                         if i.shape[0] > min_vertices]        
 
         # Output resulting lines into a single combined MultiLineString
         return MultiLineString(line_features)
@@ -537,8 +539,9 @@ def contours_to_multiline(da_i, z_value, min_vertices=2):
     if output_path and output_path.endswith('.geojson'):
         if verbose:
             print(f'Writing contours to {output_path}')
-        contours_gdf.to_crs({'init': 'EPSG:4326'}).to_file(filename=output_path, 
-                                                           driver='GeoJSON')
+        contours_gdf.to_crs('EPSG:4326').to_file(filename=output_path, 
+                                                 driver='GeoJSON')
+
     if output_path and output_path.endswith('.shp'):
         if verbose:
             print(f'Writing contours to {output_path}')