Merge recent changes to develop to stable branch (#1104)

* Update DEA_Water_Observations.ipynb (#1103)

Update doco on bit flag locations

* Feature witbook (#1102)

* WIT to feature branch

* DEA_WIT intro notebook #902

* WIT notebook pr prep

* Hird supplementary material

* DEA_WIT intro notebook #902

* small wording updates

* added updated flowchart

* update readme to add WIT notebook

* ready for PR

* checking all cells ran

* Use dask to improve load

* Update test_notebooks.sh

try to ignore the WIT notebook in testing.
This can be updated if we get a database link working in our notebook testing.

* formatted wetlands module

* normalise spelling

* includimg

* fix dyanmic sp

* Wording updates

* update comment capitalisation

* addressing comments adding TCW detail

---------

Co-authored-by: BexDunn <[email protected]>
Co-authored-by: robbibt <[email protected]>
Co-authored-by: Bex Dunn <[email protected]>

---------

Co-authored-by: James Miller <[email protected]>
Co-authored-by: LaurenSchenk1 <[email protected]>
Co-authored-by: robbibt <[email protected]>

DEA_products/DEA_Water_Observations.ipynb

@@ -830,14 +830,14 @@
     "\n",
     "| Attribute | Bit / position   | Decimal value |\n",
     "|------|------|----|\n",
-    "| No data | 0:   `0-------` or `1-------` | 1|\n",
-    "| Non contiguous | 1:   `-0------` or `-1------` | 2 |\n",
-    "| Sea | 2:   `--0-----` or `--1-----` | 4 |\n",
-    "| Terrain or low solar angle | 3:   `---0----` or `---1----` | 8 |\n",
-    "| High slope | 4:   `----0---` or `----1---` | 16 |\n",
-    "| Cloud shadow | 5:   `-----0--` or `-----1--` | 32 |\n",
-    "| Cloud | 6:   `------0-` or `------1-` | 64 |\n",
-    "| Water | 7:   `-------0` or `-------1` | 128 |\n",
+    "| No data | 0:   `-------0` or `-------1` | 1|\n",
+    "| Non contiguous | 1:   `------0-` or `------1-` | 2 |\n",
+    "| Low solar angle | 2:   `-----0--` or `-----1--` | 4 |\n",
+    "| Terrain shadow | 3:   `----0---` or `----1---` | 8 |\n",
+    "| High slope | 4:   `---0----` or `---1----` | 16 |\n",
+    "| Cloud shadow | 5:   `--0-----` or `--1-----` | 32 |\n",
+    "| Cloud | 6:   `-0------` or `-1------` | 64 |\n",
+    "| Water observed | 7:   `0-------` or `1-------` | 128 |\n",
     "\n",
     "Any combinations of flags can be combined to create a unique decimal value. \n",
     "For example, a value of 136 indicates that water (128) AND terrain shadow / low solar angle (8) were observed for the pixel (i.e. 128 + 8 = 136).\n",

DEA_products/README.rst

@@ -15,6 +15,7 @@ Notebooks introducing DEA's satellite datasets and derived products, including h
    DEA_High_and_Low_Tide_Imagery.ipynb
    DEA_Coastlines.ipynb
    DEA_Waterbodies.ipynb
+   DEA_Wetlands_Insight_Tool.ipynb
 
 Citing DEA Notebooks
 --------------------

Supplementary_data/DEA_Wetlands_Insight_Tool/hird.cpg

@@ -0,0 +1 @@
+ISO-8859-1

Supplementary_data/DEA_Wetlands_Insight_Tool/hird.prj

@@ -0,0 +1 @@
+PROJCS["GDA_1994_Australia_Albers",GEOGCS["GCS_GDA_1994",DATUM["D_GDA_1994",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Albers"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",132.0],PARAMETER["Standard_Parallel_1",-18.0],PARAMETER["Standard_Parallel_2",-36.0],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

Tests/test_notebooks.sh

@@ -11,6 +11,6 @@ pip3 install ./Tools
 pytest Tests/dea_tools
 
 # Test Juputer Notebooks
-pytest --durations=10 --nbval-lax Beginners_guide DEA_products How_to_guides/Animated_timeseries.ipynb How_to_guides/Contour_extraction.ipynb How_to_guides/Calculating_band_indices.ipynb How_to_guides/Downloading_data_with_STAC.ipynb How_to_guides/Exporting_GeoTIFFs.ipynb How_to_guides/Generating_composites.ipynb How_to_guides/Image_segmentation.ipynb How_to_guides/Opening_GeoTIFFs_NetCDFs.ipynb How_to_guides/Pansharpening.ipynb How_to_guides/Planetary_computer.ipynb How_to_guides/Polygon_drill.ipynb How_to_guides/Principal_component_analysis.ipynb How_to_guides/Rasterize_vectorize.ipynb How_to_guides/Tidal_modelling.ipynb How_to_guides/Using_load_ard.ipynb How_to_guides/Virtual_products.ipynb Real_world_examples/Coastal_erosion.ipynb Real_world_examples/Intertidal_elevation.ipynb
+pytest --durations=10 --nbval-lax Beginners_guide DEA_products --ignore DEA_products/DEA_Wetlands_Insight_Tool.ipynb How_to_guides/Animated_timeseries.ipynb How_to_guides/Contour_extraction.ipynb How_to_guides/Calculating_band_indices.ipynb How_to_guides/Downloading_data_with_STAC.ipynb How_to_guides/Exporting_GeoTIFFs.ipynb How_to_guides/Generating_composites.ipynb How_to_guides/Image_segmentation.ipynb How_to_guides/Opening_GeoTIFFs_NetCDFs.ipynb How_to_guides/Pansharpening.ipynb How_to_guides/Planetary_computer.ipynb How_to_guides/Polygon_drill.ipynb How_to_guides/Principal_component_analysis.ipynb How_to_guides/Rasterize_vectorize.ipynb How_to_guides/Tidal_modelling.ipynb How_to_guides/Using_load_ard.ipynb How_to_guides/Virtual_products.ipynb Real_world_examples/Coastal_erosion.ipynb Real_world_examples/Intertidal_elevation.ipynb
 
 

Tools/dea_tools/wetlands.py

@@ -0,0 +1,311 @@
+# wetlands.py
+"""
+This module is for processing DEA wetlands data. 
+
+License: The code in this notebook is licensed under the Apache 
+License,Version 2.0 (https://www.apache.org/licenses/LICENSE-2.0). 
+Digital Earth Australia data is licensed under the Creative Commons 
+by Attribution 4.0 license 
+(https://creativecommons.org/licenses/by/4.0/).
+
+Contact: If you need assistance, please post a question on the Open 
+Data Cube Slack channel (http://slack.opendatacube.org/) or on the 
+GIS Stack Exchange 
+(https://gis.stackexchange.com/questions/ask?tags=open-data-cube)using
+the `open-data-cube` tag (you can view previously asked questions
+here: https://gis.stackexchange.com/questions/tagged/open-data-cube). 
+
+If you would like to report an issue with this script, file one on 
+Github: https://github.com/GeoscienceAustralia/dea-notebooks/issues/new
+
+Last modified: July 2023
+"""
+
+import seaborn as sns
+import datetime
+import matplotlib.dates as mdates
+import matplotlib.pyplot as plt
+import pandas as pd
+
+# disable DeprecationWarning for chained assignments in conversion to
+# datetime format
+pd.options.mode.chained_assignment = None  # default='warn'
+
+
+def normalise_wit(polygon_base_df):
+    """
+    This function is to normalise the Fractional Cover vegetation
+    components so users can choose to display the WIT plot in a more
+    readable way. Normalising vegetation components so they total to 1.
+    Normalised values are returned as additional columns.
+
+    Last modified: July 2023
+
+    Parameters
+    ----------
+    polygon_base_df : pandas DataFrame with columns:
+    ['date',
+    'pv',
+    'npv',
+    'bs',
+    'wet',
+    'water']
+
+    Returns
+    -------
+    polygon_base_df with columns:
+    ['index',
+     'date',
+     'pv',
+     'npv',
+     'bs',
+     'wet',
+     'water',
+     'veg_areas',
+     'overall_veg_num',
+     'norm_bs',
+     'norm_pv',
+     'norm_npv']
+
+    Example
+    --------
+
+    A polygon has 11 pixels
+
+    [cloud][water][wet][wet][wet][wet][wet][wet][wet][wet][vegetation]
+      |      |        |                                        |
+      |      |        |                                        |
+      |      |        |__> wet = 8/10 = 80%                    |__> pv/npv/bs == 1/10 = 10%
+      |      |
+      |      |__> water = 1/10 = 10%
+      |
+      |__> pc_missing = 1/11 ~+ 9.1%
+
+    The vegetation pixel relative np, npv, and bs values
+
+     [vegetation]
+          |
+          |__> [pv] [npv] [bs]
+               [ 5] [  4] [ 2]
+
+    Assume vegetation relative values are:
+
+    water = 0.1
+    wet = 0.8
+
+    pv = 0.05
+    npv = 0.04
+    bs = 0.02
+
+    vegetation_area = 1 - water - wet
+
+    vegetation_overall_value = pv + npv + bs
+
+    print(f"The pv is {pv} \nThe npv is {npv} \nThe bs is {bs}
+    \nThe overall number is {water + wet + pv + npv + bs}")
+
+    The pv is 0.05
+    The npv is 0.04
+    The bs is 0.02
+    The overall number is 1.01
+
+    The overall number is greater than 1. Let us normalise the result.
+    The water and wet are pixel classification result, so we should
+    not touch them.
+
+    pv = pv/vegetation_overall_value*vegetation_area
+    npv = npv/vegetation_overall_value*vegetation_area
+    bs = bs/vegetation_overall_value*vegetation_area
+
+    print(f"The normalised pv is {pv} \nThe normalised npv is {npv}
+    \nThe normalised bs is {bs} \nThe normalised overall number is
+    {water + wet + pv + npv + bs}")
+
+    The normalised pv is 0.04545454545454545
+    The normalised npv is 0.036363636363636355
+    The normalised bs is 0.018181818181818177
+    The normalised overall number is 1.0
+
+    """
+
+    # ignore high pixel missing timestamp result
+    polygon_base_df = polygon_base_df.dropna(subset=["bs"])
+
+    # 1. compute the expected vegetation area total size: 1 - water (%) - wet (%)
+    polygon_base_df.loc[:, "veg_areas"] = (
+        1 - polygon_base_df["water"] - polygon_base_df["wet"]
+    )
+
+    # 2. normalise the vegetation values based on vegetation size (to handle FC values more than 100 issue)
+    # WARNNING: Not touch the water and wet, cause they are pixel classification result
+    polygon_base_df.loc[:, "overall_veg_num"] = (
+        polygon_base_df["pv"] + polygon_base_df["npv"] + polygon_base_df["bs"]
+    )
+
+    # 3. if the overall_veg_num is 0, no need to normalize veg area
+    norm_veg_index = polygon_base_df["overall_veg_num"] != 0
+
+    for band in ["bs", "pv", "npv"]:
+        polygon_base_df.loc[:, "norm_" + band] = polygon_base_df.loc[:, band]
+        polygon_base_df.loc[norm_veg_index, "norm_" + band] = (
+            polygon_base_df.loc[norm_veg_index, band]
+            / polygon_base_df.loc[norm_veg_index, "overall_veg_num"]
+            * polygon_base_df.loc[norm_veg_index, "veg_areas"]
+        )
+
+    # convert the string to Python datetime format, easy to do display the result in PNG
+    polygon_base_df.loc[:, "date"] = pd.to_datetime(
+        polygon_base_df["date"], infer_datetime_format=True
+    )
+
+    polygon_base_df.reset_index(inplace=True)
+
+    return polygon_base_df
+
+
+def generate_low_quality_data_periods(df):
+    """
+    This function generates low quality data periods, including the SLC off period: https://www.usgs.gov/faqs/what-landsat-7-etm-slc-data
+    and periods with an observation density of less than four observations within a twelve month (365 days) period.
+    Off value is 100 where there is low data quality and 0 for good data.
+
+    Last modified: July 2023
+
+    Parameters
+    ----------
+    df : pandas DataFrame with columns including:
+    ['date']
+
+    Returns
+    -------
+    df : pandas DataFrame with additional column:
+    ['off_value']
+
+    """
+
+    # default: all data points are good
+    df.loc[:, "off_value"] = 0
+
+    # Add the first no-data times (SLC-off only)
+    LS5_8_gap_start = datetime.datetime(2011, 11, 1)
+    LS5_8_gap_end = datetime.datetime(2013, 4, 1)
+
+    df.loc[
+        df[(df["date"] >= LS5_8_gap_start) & (df["date"] <= LS5_8_gap_end)].index,
+        "off_value",
+    ] = 100
+
+    # periods with an observation density of less than four observations within a twelve month (365 days) period
+    for i in range(3, len(df) - 3):
+        # can change to another threshold (like: 100 days) to test dynamic no-data-period display
+        if ((df.loc[i + 3, "date"] - df.loc[i, "date"]).days) > 365:
+            df.loc[
+                df[
+                    (df["date"] >= df.loc[i, "date"])
+                    & (df["date"] <= df.loc[i + 3, "date"])
+                ].index,
+                "off_value",
+            ] = 100
+
+    return df
+
+
+def display_wit_stack_with_df(
+    polygon_base_df,
+    polygon_name="your_wetland_name",
+    png_name="your_file_name",
+    width=32,
+    height=6,
+):
+    """
+    This functions produces WIT plots. Function displays a stack plot and saves as a png.
+
+    Last modified: July 2023
+
+    Parameters
+    ----------
+    polygon_base_df : pandas DataFrame with columns including:
+    ['date',
+     'wet',
+     'water',
+     'norm_bs',
+     'norm_pv',
+     'norm_npv']
+     polygon_name : string
+     png_name : string
+
+    """
+
+    plt.rcParams["axes.facecolor"] = "white"
+    plt.rcParams["savefig.facecolor"] = "white"
+    plt.rcParams["text.usetex"] = False
+
+    fig = plt.figure()
+    fig.set_size_inches(width, height)
+    ax = fig.add_subplot(111)
+    ax.autoscale(enable=True)
+
+    pal = [
+        sns.xkcd_rgb["cobalt blue"],
+        sns.xkcd_rgb["neon blue"],
+        sns.xkcd_rgb["grass"],
+        sns.xkcd_rgb["beige"],
+        sns.xkcd_rgb["brown"],
+    ]
+
+    plt.title(
+        f"Percentage of area dominated by WOfS, Wetness, Fractional Cover for\n {polygon_name}",
+        fontsize=16,
+    )
+
+    ax.stackplot(
+        polygon_base_df["date"],
+        polygon_base_df["water"] * 100,
+        polygon_base_df["wet"] * 100,
+        polygon_base_df["norm_pv"] * 100,
+        polygon_base_df["norm_npv"] * 100,
+        polygon_base_df["norm_bs"] * 100,
+        colors=pal,
+        alpha=0.7,
+    )
+
+    # manually change the legend display order
+    legend = ax.legend(
+        ["open water", "wet", "green veg", "dry veg", "bare soil"][::-1],
+        loc="lower left",
+    )
+    handles = legend.legendHandles
+
+    for i, handle in enumerate(handles):
+        handle.set_facecolor(pal[::-1][i])
+        handle.set_alpha(0.7)
+
+    # setup the display ranges
+    ax.set_ylim(0, 100)
+    ax.set_xlim(polygon_base_df["date"].min(), polygon_base_df["date"].max())
+
+    # add a new column: 'off_value' based on low quality data setting
+    polygon_base_df = generate_low_quality_data_periods(polygon_base_df)
+
+    ax.fill_between(
+        polygon_base_df["date"],
+        0,
+        100,
+        where=polygon_base_df["off_value"] == 100,
+        color="white",
+        alpha=0.5,
+        hatch="//",
+    )
+
+    # modify the xaxis settings
+    ax.xaxis.set_major_locator(mdates.YearLocator())
+    ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y"))
+
+    x_label_text = "The Fractional Cover algorithm developed by the Joint Remote Sensing Research Program and\n the Water Observations from Space algorithm developed by Geoscience Australia are used in the production of this data"
+
+    ax.set_xlabel(x_label_text, style="italic")
+
+    plt.savefig(f"{png_name}.png", bbox_inches="tight")
+    plt.show()
+
+    plt.close(fig)