<jupyter_start><jupyter_text>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  ----  SLIPSLIP 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](slip_import)* [Choose Platform and Product](slip_plat_prod)* [Define the Extents of the Analysis](slip_define_extents)* [Load Data from the Data Cube](slip_load_data)* [Change Detection](slip_change_detect)* [NDWI (Nomalized Difference Water Index)](slip_ndwi)* [RED Reflectance](slip_red)* [ASTER Global Elevation Models](slip_aster)* [Reviewing the Evolution of the SLIP Product](slip_evo)* [Visual Comparison of SLIP Output and Baseline Composited Scene](slip_compare_output_baseline)  Import Dependencies and Connect to the Data Cube [&9652;](slip_top)<jupyter_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()<jupyter_output><empty_output><jupyter_text>Choose Platform and Product [&9652;](slip_top)<jupyter_code>platform = 'LANDSAT_8'
product = 'ls8_usgs_sr_scene'
collection = 'c1'
level = 'l2'<jupyter_output><empty_output><jupyter_text>Define the Extents of the Analysis [&9652;](slip_top)<jupyter_code># 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)<jupyter_output><empty_output><jupyter_text>Load Data from the Data Cube [&9652;](slip_top)<jupyter_code># 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()<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>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()<jupyter_output><empty_output><jupyter_text>----   Change Detection [&9652;](slip_top)  In the context of SLIP, Change detection happens through the comparison of 'current' values against 'past' values.  Trivialized Example:  $$ \Delta Value = (Value_{new} - Value_{old})/ Value_{old} $$  It is easy to define NEW as the current value being analyzed.<jupyter_code>new = acq_to_show<jupyter_output><empty_output><jupyter_text>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.      The following figure illustrates such a compositing method:      ![img](../diagrams/slip/comp_00.png)   In the figure above, t4 values are the average of t1-t3 (assuming a window size of 3)  The code below composites with a window size of 5.<jupyter_code># Generate a moving average of n values leading up to current time.
baseline = generate_baseline(landsat_ds, composite_size = 3, mode = 'average')<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>(len(new.time), len(baseline.time))<jupyter_output><empty_output><jupyter_text>What this composite looks like<jupyter_code>display_at_time([baseline, new], time = time_to_show, width = 2, w = 12)<jupyter_output><empty_output><jupyter_text>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. ----     NDWI (Nomalized Difference Water Index) [&9652;](slip_top)SLIP makes the major assumption that landslides will strip a hill/mountain-side of all of its vegetation.    SLIP uses [NDWI](https://en.wikipedia.org/wiki/Normalized_difference_water_index), 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 BANDMATHNDWI 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}}$$    $$ NDWI_{BASELINE} =  \frac{NIR_{BASELINE} - SWIR_{BASELINE}}{NIR_{BASELINE} + SWIR_{BASELINE}}$$    $$\Delta NDWI = NDWI_{NEW} - NDWI_{BASELINE}$$   The code is just as simple:<jupyter_code>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<jupyter_output><empty_output><jupyter_text>Filtering NDWIIn 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%)     $$  ndwi\_filter(Dataset) = \left\{     \begin{array}{lr}       Dataset & : |  \Delta NDWI(Dataset)  | > 0.2\\       np.nan & : |  \Delta NDWI(Dataset)  | \le 0.2     \end{array}   \right.\\ $$        In code, it's even simpler:<jupyter_code>new_ndwi_filtered = new.where(abs(ndwi_change) > 0.2)<jupyter_output><empty_output><jupyter_text>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  .<jupyter_code>display_at_time([new, (new, new_ndwi_filtered),new_ndwi_filtered], 
                time = time_to_show,
                width = 3, w =14)<jupyter_output><empty_output><jupyter_text>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.     ----   RED Reflectance [&9652;](slip_top)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 bandmathRed change is computed on a per-pixel level and involves arithmetic on the RED band values.  The derivation of RED change is simple:  $$ \Delta Red = \frac{RED_{NEW} - RED_{BASELINE}}{RED_{BASELINE}} $$ The code is just as simple:<jupyter_code>red_change = (new.red - baseline.red)/(baseline.red)<jupyter_output><empty_output><jupyter_text>Filtering for RED reflectance increaseFiltering RED reflectance change is just like the piecewise transformation used for filtering NDWI change.      $$  red\_filter(Dataset) = \left\{     \begin{array}{lr}       Dataset & :   \Delta red(Dataset)   > 0.4\\       np.nan & :   \Delta red(Dataset)   \le 0.4     \end{array}   \right.\\ $$  In Code:<jupyter_code>new_red_and_ndwi_filtered = new_ndwi_filtered.where(red_change > 0.4)<jupyter_output><empty_output><jupyter_text>How much further RED reflectance filtering gets youContinuing 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.<jupyter_code>display_at_time([new, (new, new_red_and_ndwi_filtered),new_red_and_ndwi_filtered], 
                time = time_to_show, 
                width = 3, w = 14)<jupyter_output><empty_output><jupyter_text>---- ASTER Global Elevation Models [&9652;](slip_top)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<jupyter_code>aster = dc.load(product="terra_aster_gdm",\
                lat=lat,\
                lon=lon,\
                measurements=['dem'],
                group_by='solar_day')<jupyter_output><empty_output><jupyter_text>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. $$ Gradient = \frac{Rise}{Run}  $$  Basic trigonometric identities  can then be used to derive the angle:       $$ Angle of Elevation = \arctan(Gradient) $$        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                ![img](diagrams/slip/derivation.png)        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.<jupyter_code># 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)<jupyter_output><empty_output><jupyter_text>Filtering out pixels that don't meet requirements for steepnessFiltering based on slope is a peicewise transformation using a derived slopemask:         $$  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.\\ $$        Its use in code:<jupyter_code>slip_product = new_red_and_ndwi_filtered.where(is_above_slope_threshold)<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>display_at_time([new, (new, slip_product),slip_product], 
                time = time_to_show,
                width = 3, w = 14)<jupyter_output><empty_output><jupyter_text>Reviewing the Evolution of the SLIP Product [&9652;](slip_top)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<jupyter_code>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)<jupyter_output><empty_output><jupyter_text>Visual Comparison of SLIP Output and Baseline Composited Scene [&9652;](slip_top)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.<jupyter_code>display_at_time([baseline, (new,slip_product)], 
                time = time_to_show, 
                width = 2, mode = 'blend', color = [210,7,7] , w = 14)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08478260869565217, 'fraction_numerical': 0.012111801242236025, 'mean_word_length': 4.534891715366105, 'pattern_counts': {'":': 0, '<': 86, '<?xml version=': 0, '>': 91, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '315611', 'license': 'Apache-2.0', 'n_tokens_mistral': 5045, 'n_tokens_neox': 4738, 'n_words': 1703, 'path': 'notebooks/landslides/Landslide_Identification_SLIP.ipynb', 'repo_name': 'admariner/data_cube_notebooks', 'size': '16100'}