In [1]:
from sys import path
path.append("../")
%matplotlib inline  
import matplotlib.pyplot as plt

Coastline Classifier

This coastal boundary algorithm is used to classify a given pixel as either coastline or not coastline using a simple binary format like in the table before.


$\begin{array}{|c|c|} \hline 1& Coastline \\ \hline 0& Not Coastline \\ \hline \end{array}$


The algorithm makes a classification by examining surrounding pixels and making a determination based on how many pixels around it are water



If the count of land pixels surrounding a pixel exceeds 5, then it's likely not coastline. If the count of land pixels surrounding a pixel does not exceed 1, then it's likely not a coastline


$$ Classification(pixel) = \begin{cases} 1 & 2\le count\_water\_surrounding(pixel) \leq 5 \\ 0 & \end{cases} $$


Counting by applying a convolutional kernel

A convolution applies a kernel to a point and it's surrounding pixels. Then maps the product to a new grid.

In the case of coastal boundary classification, A convolution the following kernel is applied to a grid of water, not-water pixels.


$$ Kernel = \begin{bmatrix} 1 & 1 & 1\\ 1 & 0 & 1\\ 1 & 1 & 1\\ \end{bmatrix} $$



There exist more complicated differential kernels that would also work( see sobel operator).
The one used in this notebooks operates on binary variables and is easier to work with and easy to debug.


Index

Import Dependencies and Connect to the Data Cube [▴](#coastline_classifier_top)


In [2]:
import scipy.ndimage.filters as conv
import numpy as np

def _coastline_classification(dataset, water_band='wofs'):
    kern = np.array([[1, 1, 1], [1, 0.001, 1], [1, 1, 1]])
    convolved = conv.convolve(dataset[water_band], kern, mode='constant') // 1

    ds = dataset.where(convolved > 0)
    ds = ds.where(convolved < 6)
    ds.wofs.values[~np.isnan(ds.wofs.values)] = 1
    ds.wofs.values[np.isnan(ds.wofs.values)] = 0

    return ds.rename({"wofs": "coastline"})

In [3]:
import datacube
dc = datacube.Datacube(app = "Coastline classification")

Choose Platform and Product [▴](#coastline_classifier_top)


In [4]:
platform        = 'LANDSAT_8'
product_type    = 'ls8_usgs_sr_scene'

Define the Extents of the Analysis [▴](#coastline_classifier_top)

West Africa is subject to considerable coastal erosion in some areas. The links listed below are references regarding coastal erosion in West Africa and coastal erosion in general.

  • World Bank WACA program brochure (2015) - link
  • USAID - Adapting to Coastal Climate Change (2009) - - link

In [5]:
# Ghana
lon = (0.0520, 0.3458)
lat = (5.6581, 5.8113)

Visualize the selected area


In [6]:
from utils.data_cube_utilities.dc_display_map import display_map
display_map(lat, lon)


Out[6]:

Load Data from the Data Cube and Create a Composite [▴](#coastline_classifier_top)


In [7]:
from datetime import datetime 

params = dict(platform=platform,
              product=product_type,
              time=(datetime(2013,1,1), datetime(2013,12,31)) ,
              lon= lon,
              lat= lat,
              measurements = ['red', 'green', 'blue', 'nir', 'swir1', 'swir2', 'pixel_qa'] )

In [8]:
dataset = dc.load(**params)

Obtain the clean mask


In [9]:
from utils.data_cube_utilities.clean_mask import landsat_qa_clean_mask, landsat_clean_mask_invalid

clean_mask = (landsat_qa_clean_mask(dataset, platform) &
             ((dataset != -9999).to_array().all('variable')) &
              landsat_clean_mask_invalid(dataset)).persist()

Create a composite


In [10]:
# from utils.data_cube_utilities.clean_mask import landsat_qa_clean_mask
# from utils.data_cube_utilities.dc_mosaic import create_median_mosaic

# def mosaic(dataset):
#     water_and_land_mask = landsat_qa_clean_mask(dataset, platform, 
#                                                 cover_types=['clear', 'water'])
#     return create_median_mosaic(dataset, clean_mask = water_and_land_mask)

In [11]:
# np.unique(clean_mask.red, return_counts=True)

In [12]:
from utils.data_cube_utilities.dc_mosaic import create_median_mosaic
from utils.data_cube_utilities.dc_utilities import ignore_warnings

composited_dataset = ignore_warnings(create_median_mosaic, dataset, clean_mask)

Visualize Composited imagery


In [13]:
from utils.data_cube_utilities.plotter_utils import figure_ratio

composited_dataset.swir1.plot(cmap = "Greys", figsize = figure_ratio(dataset, fixed_width = 20))


Out[13]:
<matplotlib.collections.QuadMesh at 0x7f514ad073c8>

Obtain Water Classifications and Coastal Change [▴](#coastline_classifier_top)


In [14]:
from utils.data_cube_utilities.dc_water_classifier import wofs_classify

water_classification = ignore_warnings(wofs_classify, composited_dataset, mosaic = True)

In [15]:
water_classification.wofs.plot(cmap = "Blues", figsize = figure_ratio(dataset, fixed_width = 20))


Out[15]:
<matplotlib.collections.QuadMesh at 0x7f516ed08358>



In [16]:
coast = _coastline_classification(water_classification, water_band='wofs')

In [17]:
coast.coastline.plot(cmap = "Blues", figsize = figure_ratio(dataset, fixed_width = 20))


Out[17]:
<matplotlib.collections.QuadMesh at 0x7f516ecbc160>