<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Sudden Landslide Identification Product (SLIP)
What to expect from this notebook
Introduction to the SLIP algorithm 
describing change detection in the context of datacube
Detailed band math equations for SLIP filtering  
Illustrate the step by step evolution of a SLIP product  
<a id='slip_top'></a>
SLIP
SLIP is used to automate the detection of Landslides. A SLIP product is the result of filtering based on per-pixel changes in both soil moisture and vegetation in areas with high elevation gradients.  All of which (with the exception of elevation gradients) can be computed using simple bandmath equations.
Data
SLIP makes use of the following Landsat 7 Surface Reflectance Bands
Step1:  <span id="slip_plat_prod">Choose Platform and Product &#9652;</span>
Step2:  <span id="slip_define_extents">Define the Extents of the Analysis &#9652;</span>
Step3:  <span id="slip_load_data">Load Data from the Data Cube &#9652;</span>
Step4:  Visualization
This step is optional, but useful to those seeking a step by step validation of SLIP. The following code shows a true-color representation of our loaded scene.
Step5:  <span id="slip_change_detect">Change Detection &#9652;</span>
In the context of SLIP, Change detection happens through the comparison of 'current' values against 'past' values.
<br>
Trivialized Example
Step6:  <br>
However, OLD can have varying interpretations.
In SLIP, OLD values (referred to in code as BASELINE values) are simply rolling averages of not-nan values leading up to the date in question.  
<br>
The following figure illustrates such a compositing method
Step7:  It is important to note that compositing will shorten the length of baseline's time domain by the window size since ranges less than the composite size are not computed. For a composite size of 5, new's first 5 time values will not have composite values.
Step8:  What this composite looks like
Step9:  The baseline composite is featured in the figure above (left). It represents what was typical for the past five acquisitions 'leading-up-to' time_to_show. Displayed next to it (right) is the true-color visualization of the acquisition 'at' time_to_show. The new object contains unaltered LS7 scenes that are index-able  using a date like time_to_show. The baseline object contains a block of composites of those landsat scenes that is index-able the same way.
<span id="slip_ndwi">NDWI (Nomalized Difference Water Index) &#9652;</span>
SLIP makes the major assumption that landslides will strip a hill/mountain-side of all of its vegetation.  
SLIP uses NDWI, an index used to monitor water content of leaves, to track the existence of vegetation on a slope. At high enough levels, leaf water content change can no longer be attributed to something like seasonal fluctuations and will most likely indicate a change in the existence of vegetation.  
NDWI BANDMATH
NDWI is computed on a per-pixel level and involves arithmetic between NIR (Near infrared) and SWIR1 (Short Wave Infrared) values. 
NDWI is computed for both NEW and BASELINE imagery then compared to yield NDWI change. The equations bellow detail a very simple derivation of change in NDWI
Step10:  Filtering NDWI
In the context of code, you can best think of filtering as a peicewise transformation that assigns a nan (or null) value to points that fall below our minimum change threshold. (For SLIP that threshold is 20%) 
<br>
$$  ndwi_filter(Dataset) = \left{
     \begin{array}{lr}
       Dataset & 
Step11:  How far NDWI filtering gets you
A SLIP product is the result of a process of elimination.  NDWI is sufficient in eliminating a majority of non-contending areas early on in the process. Featured below is what is left of the original image after having filtered for changes in NDWI  .
Step12:  Highlighted in the center picture are values that meet our NDWI change expectations. Featured in the right-most image is what remains of our original image after NDWI filtering.  
<span id="slip_red">RED Reflectance &#9652;</span>
SLIP makes another important assumption about Landslides. 
On top of stripping the Slope of vegetation, a landslide will reveal a large layer of previously vegetated soil. Since soil reflects more light in the RED spectral band than highly vegetated areas do, SLIP looks for increases in the RED bands. This captures both the loss of vegetation, and the unearthing of soil.    
RED change bandmath
Red change is computed on a per-pixel level and involves arithmetic on the RED band values.  The derivation of RED change is simple
Step13:  Filtering for RED reflectance increase
Filtering RED reflectance change is just like the piecewise transformation used for filtering NDWI change.  
<br>
$$  red_filter(Dataset) = \left{
     \begin{array}{lr}
       Dataset & 
Step14:  How much further RED reflectance filtering gets you
Continuing SLIP's process of elimination, Red increase filtering will further refine the area of interest to areas that, upon visual inspection appear to be light brown in color.
Step15:  <span id="slip_aster">ASTER Global Elevation Models &#9652;</span>
Aster GDEM models provide elevation data for each pixel expressed in meters.  For SLIP height is not enough to determine that a landslide can happen on a pixel. SLIP focuses on areas with high elevation Gradients/Slope (Expressed in non-radian degrees).The driving motivation for using slope based filtering is that landslides are less likely to happen in flat regions. 
Loading the elevation model
Step16:  Calculating Angle of elevation
A gradient is generated for each pixel using the four pixels adjacent to it, as well as a rise/run formuala. 
<br><br>
$$ Gradient = \frac{Rise}{Run}  $$
<br><br>
Basic trigonometric identities  can then be used to derive the angle
Step17:  Filtering out pixels that don't meet requirements for steepness
<br>
Filtering based on slope is a peicewise transformation using a derived slopemask
Step18:  Visualising our final SLIP product
The final results of SLIP are small regions of points with a high likelihood of landslides having occurred on them. Furthermore there is no possibility that detections are made in flat areas(areas with less than a $15^{\circ}$ angle of elevation.
Step19:  <span id="slip_evo">Reviewing the Evolution of the SLIP Product &#9652;</span>
The following visualizations will detail the evolution of the SLIP product from the previous steps.
Order of operations
Step20:  <span id="slip_compare_output_baseline">Visual Comparison of SLIP Output and Baseline Composited Scene &#9652;</span>
In the name of validating results, it makes sense to compare the SLIP product generated for the selected date (time_to_show) to the composited scene representing what is considered to be "normal" for the last 5 acquisitions.
<ASSISTANT_TASK:>
Python Code:
import sys
import os
sys.path.append(os.environ.get('NOTEBOOK_ROOT'))
import numpy as np
import xarray as xr
import pandas as pd
import matplotlib.pyplot as plt
from utils.data_cube_utilities.dc_display_map import display_map
from utils.data_cube_utilities.clean_mask import landsat_clean_mask_full
#     landsat_qa_clean_mask, landsat_clean_mask_invalid
from utils.data_cube_utilities.dc_baseline import generate_baseline
from utils.data_cube_utilities.dc_displayutil import display_at_time
from utils.data_cube_utilities.dc_slip import create_slope_mask
from datacube.utils.aws import configure_s3_access
configure_s3_access(requester_pays=True)
import datacube  
dc = datacube.Datacube()
Explanation: Sudden Landslide Identification Product (SLIP)
What to expect from this notebook
Introduction to the SLIP algorithm 
describing change detection in the context of datacube
Detailed band math equations for SLIP filtering  
Illustrate the step by step evolution of a SLIP product  
<a id='slip_top'></a>
SLIP
SLIP is used to automate the detection of Landslides. A SLIP product is the result of filtering based on per-pixel changes in both soil moisture and vegetation in areas with high elevation gradients.  All of which (with the exception of elevation gradients) can be computed using simple bandmath equations.
Data
SLIP makes use of the following Landsat 7 Surface Reflectance Bands:
- RED,
- NIR,
- SWIR1
- PIXEL_QA  
SLIP makes use of the following ASTER GDEM V2 bands:
- dem  
Algorithmic Process
Algorithmically speaking, SLIP is a series of per-pixel filter operations acting on relationships between NEW(current) and BASELINE(historical) values of an area. The remaining pixels after filter operations will be what SLIP classifies as landslides.  Itemized in the list below are operations taken to create a SLIP product:  
Import and initialize datacube
Load Geographic area
Remove clouds and no-data values
Label this product NEW
Generate a rolling average composite of NEW  
Label the rolling average composite BASELINE
Filter in favor of sufficiently large changes in vegetation (using NDWI values derived from NEW and BASELINE)  
Filter in favor of sufficiently large increases in RED reflectance(using RED band values from NEW and BASELINE)  
Generate a slope-mask(using ASTERDEM V2 data) 
Filter in favor of areas that have a high enough slope(Landslides don't happen on flat surfaces)   
Index
Import Dependencies and Connect to the Data Cube
Choose Platform and Product
Define the Extents of the Analysis
Load Data from the Data Cube
Change Detection
NDWI (Nomalized Difference Water Index)
RED Reflectance
ASTER Global Elevation Models
Reviewing the Evolution of the SLIP Product
Visual Comparison of SLIP Output and Baseline Composited Scene
<span id="slip_import">Import Dependencies and Connect to the Data Cube &#9652;</span>
End of explanation
platform = 'LANDSAT_8'
product = 'ls8_usgs_sr_scene'
collection = 'c1'
level = 'l2'
Explanation: <span id="slip_plat_prod">Choose Platform and Product &#9652;</span>
End of explanation
# Freetown, Sierra Leone
# (https://www.reuters.com/article/us-leone-mudslide-africa/cities-across-africa-face-threat-of-landslides-like-sierra-leone-idUSKCN1AY115)
# define geographic boundaries in (min, max) format
lon = (-13.3196, -12.9366)
lat = (8.1121, 8.5194)
# define date range boundaries in (min,max) format
# There should be a landslide by Freetown during August 2017.
date_range =("2016-01-01", "2017-12-31")
display_map(lat, lon)
Explanation: <span id="slip_define_extents">Define the Extents of the Analysis &#9652;</span>
End of explanation
# Define desired bands. For SLIP, only red, nir, swir and pixel_qa will be necessary.
desired_bands = ['red','nir','swir1','pixel_qa']  
  
# Add blue and green bands since they are needed for visualizing results (RGB).
desired_bands = desired_bands + ['green', 'blue']
  
# Load area.
landsat_ds = dc.load(product = product,\
                     platform = platform,\
                     lat = lat,\
                     lon = lon,\
                     time = date_range,\
                     measurements = desired_bands,
                     group_by='solar_day',
                     dask_chunks={'time':1, 'longitude': 1000, 'latitude': 1000}).persist()
# clean_mask = landsat_qa_clean_mask(landsat_ds, platform) & \
#              (landsat_ds != -9999).to_array().all('variable') & \
#              landsat_clean_mask_invalid(landsat_ds)
clean_mask = landsat_clean_mask_full(dc, landsat_ds, product=product, platform=platform,
                                     collection=collection, level=level).persist()
# Determine the times with data.
data_time_mask = (clean_mask.sum(['latitude', 'longitude']) > 0).persist()
clean_mask = clean_mask.sel(time=data_time_mask)
landsat_ds = landsat_ds.sel(time=data_time_mask)
landsat_ds = landsat_ds.where(clean_mask).persist()
Explanation: <span id="slip_load_data">Load Data from the Data Cube &#9652;</span>
End of explanation
time_to_show = '2017-08-04'
acq_to_show = landsat_ds.sel(time=time_to_show, method='nearest')
rgb_da = acq_to_show[['red', 'green', 'blue']].squeeze().to_array().compute()
vmin = rgb_da.quantile(0.05).values
vmax = rgb_da.quantile(0.95).values
rgb_da.plot.imshow(vmin=vmin, vmax=vmax)
plt.show()
Explanation: Visualization
This step is optional, but useful to those seeking a step by step validation of SLIP. The following code shows a true-color representation of our loaded scene.
End of explanation
new = acq_to_show
Explanation: <span id="slip_change_detect">Change Detection &#9652;</span>
In the context of SLIP, Change detection happens through the comparison of 'current' values against 'past' values.
<br>
Trivialized Example:
<br>
$$ \Delta Value = (Value_{new} - Value_{old})/ Value_{old} $$  
<br>
It is easy to define NEW as the current value being analyzed. 
<br>
End of explanation
# Generate a moving average of n values leading up to current time.
baseline = generate_baseline(landsat_ds, composite_size = 3, mode = 'average')  
Explanation: <br>
However, OLD can have varying interpretations.
In SLIP, OLD values (referred to in code as BASELINE values) are simply rolling averages of not-nan values leading up to the date in question.  
<br>
The following figure illustrates such a compositing method:  
<br><br>  
<!-- ![img](avg_compositing.png) -->
<br>
In the figure above, t4 values are the average of t1-t3 (assuming a window size of 3)  
<br>
The code below composites with a window size of 5.
End of explanation
(len(new.time), len(baseline.time))
Explanation: It is important to note that compositing will shorten the length of baseline's time domain by the window size since ranges less than the composite size are not computed. For a composite size of 5, new's first 5 time values will not have composite values.
End of explanation
display_at_time([baseline, new], time = time_to_show, width = 2, w = 12)
Explanation: What this composite looks like
End of explanation
ndwi_new = (new.nir- new.swir1)/(new.nir + new.swir1)
ndwi_baseline = (baseline.nir - baseline.swir1)/ (baseline.nir + baseline.swir1)
ndwi_change = ndwi_new - ndwi_baseline
Explanation: The baseline composite is featured in the figure above (left). It represents what was typical for the past five acquisitions 'leading-up-to' time_to_show. Displayed next to it (right) is the true-color visualization of the acquisition 'at' time_to_show. The new object contains unaltered LS7 scenes that are index-able  using a date like time_to_show. The baseline object contains a block of composites of those landsat scenes that is index-able the same way.
<span id="slip_ndwi">NDWI (Nomalized Difference Water Index) &#9652;</span>
SLIP makes the major assumption that landslides will strip a hill/mountain-side of all of its vegetation.  
SLIP uses NDWI, an index used to monitor water content of leaves, to track the existence of vegetation on a slope. At high enough levels, leaf water content change can no longer be attributed to something like seasonal fluctuations and will most likely indicate a change in the existence of vegetation.  
NDWI BANDMATH
NDWI is computed on a per-pixel level and involves arithmetic between NIR (Near infrared) and SWIR1 (Short Wave Infrared) values. 
NDWI is computed for both NEW and BASELINE imagery then compared to yield NDWI change. The equations bellow detail a very simple derivation of change in NDWI: 
$$ NDWI_{NEW} =  \frac{NIR_{NEW} - SWIR_{NEW}}{NIR_{NEW} + SWIR_{NEW}}$$
  <br><br>
$$ NDWI_{BASELINE} =  \frac{NIR_{BASELINE} - SWIR_{BASELINE}}{NIR_{BASELINE} + SWIR_{BASELINE}}$$
  <br><br>
$$\Delta NDWI = NDWI_{NEW} - NDWI_{BASELINE}$$  
<br>
The code is just as simple:
End of explanation
new_ndwi_filtered = new.where(abs(ndwi_change) > 0.2)
Explanation: Filtering NDWI
In the context of code, you can best think of filtering as a peicewise transformation that assigns a nan (or null) value to points that fall below our minimum change threshold. (For SLIP that threshold is 20%) 
<br>
$$  ndwi_filter(Dataset) = \left{
     \begin{array}{lr}
       Dataset & : |  \Delta NDWI(Dataset)  | > 0.2\
       np.nan & : |  \Delta NDWI(Dataset)  | \le 0.2
     \end{array}
   \right.\ $$  
<br>
In code, it's even simpler:
End of explanation
display_at_time([new, (new, new_ndwi_filtered),new_ndwi_filtered], 
                time = time_to_show,
                width = 3, w =14)
Explanation: How far NDWI filtering gets you
A SLIP product is the result of a process of elimination.  NDWI is sufficient in eliminating a majority of non-contending areas early on in the process. Featured below is what is left of the original image after having filtered for changes in NDWI  .
End of explanation
red_change = (new.red - baseline.red)/(baseline.red)
Explanation: Highlighted in the center picture are values that meet our NDWI change expectations. Featured in the right-most image is what remains of our original image after NDWI filtering.  
<span id="slip_red">RED Reflectance &#9652;</span>
SLIP makes another important assumption about Landslides. 
On top of stripping the Slope of vegetation, a landslide will reveal a large layer of previously vegetated soil. Since soil reflects more light in the RED spectral band than highly vegetated areas do, SLIP looks for increases in the RED bands. This captures both the loss of vegetation, and the unearthing of soil.    
RED change bandmath
Red change is computed on a per-pixel level and involves arithmetic on the RED band values.  The derivation of RED change is simple: 
<br><br>
$$ \Delta Red = \frac{RED_{NEW} - RED_{BASELINE}}{RED_{BASELINE}} $$
The code is just as simple:
End of explanation
new_red_and_ndwi_filtered = new_ndwi_filtered.where(red_change > 0.4)
Explanation: Filtering for RED reflectance increase
Filtering RED reflectance change is just like the piecewise transformation used for filtering NDWI change.  
<br>
$$  red_filter(Dataset) = \left{
     \begin{array}{lr}
       Dataset & :   \Delta red(Dataset)   > 0.4\
       np.nan & :   \Delta red(Dataset)   \le 0.4
     \end{array}
   \right.\ $$  
<br>
In Code:
End of explanation
display_at_time([new, (new, new_red_and_ndwi_filtered),new_red_and_ndwi_filtered], 
                time = time_to_show, 
                width = 3, w = 14)
Explanation: How much further RED reflectance filtering gets you
Continuing SLIP's process of elimination, Red increase filtering will further refine the area of interest to areas that, upon visual inspection appear to be light brown in color.
End of explanation
aster = dc.load(product="terra_aster_gdm",\
                lat=lat,\
                lon=lon,\
                measurements=['dem'],
                group_by='solar_day')  
Explanation: <span id="slip_aster">ASTER Global Elevation Models &#9652;</span>
Aster GDEM models provide elevation data for each pixel expressed in meters.  For SLIP height is not enough to determine that a landslide can happen on a pixel. SLIP focuses on areas with high elevation Gradients/Slope (Expressed in non-radian degrees).The driving motivation for using slope based filtering is that landslides are less likely to happen in flat regions. 
Loading the elevation model
End of explanation
# Create a slope-mask. False: if pixel <15 degees; True: if pixel > 15 degrees;  
is_above_slope_threshold = create_slope_mask(aster, degree_threshold = 15,resolution = 30) 
Explanation: Calculating Angle of elevation
A gradient is generated for each pixel using the four pixels adjacent to it, as well as a rise/run formuala. 
<br><br>
$$ Gradient = \frac{Rise}{Run}  $$
<br><br>
Basic trigonometric identities  can then be used to derive the angle: 
<br><br>
$$ Angle of Elevation = \arctan(Gradient) $$
   <br><br>  
When deriving the angle of elevation for a pixel, two gradients are available. One formed by the bottom pixel and top pixel, the other formed by the right and left pixel. For the purposes of identifying landslide causing slopes, the greatest of the two slopes will be used.
The following image describes the process for angle-of-elevation calculation for a single pixel within a grid of DEM pixels      
<br><br> 
<br><br>  
The vagaries of implementation have been abstracted away by dc_demutils. It's used to derive a slope-mask. A slope-mask in this sense, is an array of true and false values based on whether or not that pixel meets a minimum angle of elevation requirement.  Its use is detailed below.
End of explanation
slip_product = new_red_and_ndwi_filtered.where(is_above_slope_threshold)
Explanation: Filtering out pixels that don't meet requirements for steepness
<br>
Filtering based on slope is a peicewise transformation using a derived slopemask:   
<br>  
$$  slope_filter(Dataset) = \left{
     \begin{array}{lr}
       Dataset & :   is_above_degree_threshold(Dataset, 15^{\circ})   = True\
       np.nan &  :    is_above_degree_threshold(Dataset, 15^{\circ})   = False\
     \end{array}
   \right.\ $$  
<br>
Its use in code:
End of explanation
display_at_time([new, (new, slip_product),slip_product], 
                time = time_to_show,
                width = 3, w = 14)
Explanation: Visualising our final SLIP product
The final results of SLIP are small regions of points with a high likelihood of landslides having occurred on them. Furthermore there is no possibility that detections are made in flat areas(areas with less than a $15^{\circ}$ angle of elevation.
End of explanation
display_at_time([new, (new,new_ndwi_filtered),new_ndwi_filtered,new,  (new, new_red_and_ndwi_filtered),new_red_and_ndwi_filtered, new, (new, slip_product),slip_product], 
                time = time_to_show,
                width = 3, w = 14, h = 12)
Explanation: <span id="slip_evo">Reviewing the Evolution of the SLIP Product &#9652;</span>
The following visualizations will detail the evolution of the SLIP product from the previous steps.
Order of operations:
- NDWI change Filtered
- RED increase Filtered
- Slope Filtered
Visualization
End of explanation
display_at_time([baseline, (new,slip_product)], 
                time = time_to_show, 
                width = 2, mode = 'blend', color = [210,7,7] , w = 14)
Explanation: <span id="slip_compare_output_baseline">Visual Comparison of SLIP Output and Baseline Composited Scene &#9652;</span>
In the name of validating results, it makes sense to compare the SLIP product generated for the selected date (time_to_show) to the composited scene representing what is considered to be "normal" for the last 5 acquisitions.
End of explanation
<END_TASK>