In [1]:
from sys import path
path.append("../")
%matplotlib inline
import matplotlib.pyplot as plt
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
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.
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.
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")
In [4]:
platform = 'LANDSAT_8'
product_type = 'ls8_usgs_sr_scene'
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]:
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]:
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]:
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]: