International agencies like World Bank, UNEP, and USAID are currently reporting and addressing the problem of coastal erosion near Lomé, Togo. The links listed below are references from these agencies regarding coastal erosion in Togo and coastal erosion as a world wide phenomena.
y1 and y2, y1 and y2 to yield y1_composite and y2_composite .(This removes clouds, scanline errors, and statistical anomalies) y1_composite and y2_composite, to yield y1_water and y2_water. These outputs contain a 1 for water, and 0 for not water. y1_water from y2_water to yield a coastal_change product. Where 0 represents no change, 1 represents water gain, and -1 represents water loss. Flow Diagram
In [1]:
%matplotlib inline
from datetime import datetime
import numpy as np
import utils.data_cube_utilities.dc_utilities as utils
from utils.data_cube_utilities.dc_mosaic import create_mosaic, ls8_unpack_qa
from utils.data_cube_utilities.dc_water_classifier import wofs_classify
# Initialize data cube object
import datacube
dc = datacube.Datacube(app='dc-coastal-erosion')
In [2]:
start_time = datetime.now()
print("Start time: " + str(start_time))
In [3]:
# Set query parameters
platform = 'LANDSAT_8'
product_type = 'ls8_usgs_sr_scene'
In [4]:
# Select minimum and maximum longitudes and latitudes and the time range.
# Ghana
# Coastline east of Accra
lon = (0.0520, 0.3458)
lat = (5.6581, 5.8113)
first_year = 2013
first_time_range = (f'{first_year}-01-01', f'{first_year}-12-31')
last_year = 2014
second_time_range = (f'{last_year}-01-01', f'{last_year}-12-31')
Visualize the selected area
In [5]:
from utils.data_cube_utilities.dc_display_map import display_map
display_map(lat, lon)
Out[5]:
In [6]:
measurements = ['red', 'green', 'blue', 'nir', 'swir1', 'swir2', 'pixel_qa']
# Retrieve data from Data Cube
first_dataset = dc.load(platform=platform,
product=product_type,
time=first_time_range,
lon=lon, lat=lat, measurements=measurements)
In [7]:
# Retrieve data from Data Cube
second_dataset = dc.load(platform=platform,
product=product_type,
time=second_time_range,
lon=lon, lat=lat, measurements=measurements)
In [8]:
from utils.data_cube_utilities.dc_utilities import ignore_warnings
# Only keep pixels that are clear or have water.
clear_xarray = ls8_unpack_qa(first_dataset.pixel_qa, "clear")
water_xarray = ls8_unpack_qa(first_dataset.pixel_qa, "water")
first_clean_mask = np.logical_or(clear_xarray.astype(bool),
water_xarray.astype(bool))
clear_xarray = ls8_unpack_qa(second_dataset.pixel_qa, "clear")
water_xarray = ls8_unpack_qa(second_dataset.pixel_qa, "water")
second_clean_mask = np.logical_or(clear_xarray.astype(bool),
water_xarray.astype(bool))
# Remove noise from images by using appropriate data within the dataset to replace "dirty" data.
first_mosaic = ignore_warnings(create_mosaic, first_dataset, clean_mask=first_clean_mask)
second_mosaic = ignore_warnings(create_mosaic, second_dataset, clean_mask=second_clean_mask)
In [9]:
first_water_class = \
ignore_warnings(wofs_classify, first_mosaic, mosaic=True, no_data=np.nan)
second_water_class = \
ignore_warnings(wofs_classify, second_mosaic, mosaic=True, no_data=np.nan)
first_wofs = first_water_class.wofs.values
second_wofs = second_water_class.wofs.values
In [10]:
coastal_change = second_water_class - first_water_class
In [11]:
import matplotlib.pyplot as plt
def plot_data_array_with_aspect(da):
fig, ax = plt.subplots(figsize=(15,8))
ax.set_aspect('equal')
da.plot()
Show the coastal change
In [12]:
# -1 -> water to coast
# 0 -> no change
# 1 -> coast to water (Coastal Erosion)
plot_data_array_with_aspect( coastal_change.wofs )
Obtain the coastlines for the first and second time periods
In [13]:
first_coastline = np.zeros(first_wofs.shape)
for i in range(first_wofs.shape[0]):
for j in range(first_wofs.shape[1]):
pixel = first_wofs[i,j]
if pixel == 0 and np.nansum(first_wofs[i-1:i+2, j-1:j+2]) >= 1 and np.nansum(first_wofs[i-1:i+2, j-1:j+2]) <= 5:
first_wofs[i,j] = 1
second_coastline = np.zeros(second_wofs.shape)
for i in range(second_wofs.shape[0]):
for j in range(second_wofs.shape[1]):
pixel = second_wofs[i,j]
if pixel == 0 and np.nansum(second_wofs[i-1:i+2, j-1:j+2]) >= 1 and np.nansum(second_wofs[i-1:i+2, j-1:j+2]) <= 5:
second_wofs[i,j] = 1
Show the water classifications for the first and second time periods
In [14]:
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
In [15]:
from pylab import imshow
In [16]:
fig = plt.figure(figsize = (15,8))
a=fig.add_subplot(1,2,1)
imgplot = plt.imshow(first_wofs, cmap='Blues',
extent=[first_water_class.longitude.values.min(),
first_water_class.longitude.values.max(),
first_water_class.latitude.values.min(),
first_water_class.latitude.values.max()])
a.set_title(f'wofs - {first_year}')
plt.colorbar(ticks=[0,1], orientation ='horizontal')
a=fig.add_subplot(1,2,2)
imgplot = plt.imshow(second_wofs, cmap='Blues',
extent=[second_water_class.longitude.values.min(),
second_water_class.longitude.values.max(),
second_water_class.latitude.values.min(),
second_water_class.latitude.values.max()])
#imgplot.set_clim(0.0,1.0)
a.set_title(f'wofs - {last_year}')
plt.colorbar(ticks=[0,1], orientation='horizontal')
plt.savefig('wofs_compare.png')
Show the coastline for the first time period
In [17]:
fig = plt.figure(figsize = (15,8))
plt.imshow(first_coastline, cmap='Greens')
Out[17]:
Show the coastline for the second time period
In [18]:
fig = plt.figure(figsize = (15,8))
plt.imshow(second_coastline, cmap='Purples')
Out[18]: