Bathymetry is the measurement of depth in bodies of water(Oceans, Seas or Lakes). This notebook illustrates a technique for deriving depth of shallow water areas using purely optical features from Landsat Collection 1 imagery and draws heavily from the publication Determination of water depth with high-resolution satelite imagery over variable bottom types.
Bathymetry Index
This bathymetry index uses optical green and blue values on a logarithmic scale with two tunable coefficients m0 and m1.
Where:
m0 is a tunable scaling factor to tune the ratio to depth m1 is the offset for a depth of 0 meters.
In [1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns
def bathymetry_index(df, m0 = 1, m1 = 0):
return m0*(np.log(df.blue)/np.log(df.green))+m1
In [2]:
import datacube
dc = datacube.Datacube()
In [3]:
#List the products available on this server/device
dc.list_products()
Out[3]:
In [4]:
#create a list of the desired platforms
platform = "LANDSAT_8"
product = "ls8_level1_usgs"
In [5]:
# East Coast of Australia
lat_subsect = (-31.7, -32.2)
lon_subsect = (152.4, 152.9)
In [6]:
print('''
Latitude:\t{0}\t\tRange:\t{2} degrees
Longitude:\t{1}\t\tRange:\t{3} degrees
'''.format(lat_subsect,
lon_subsect,
max(lat_subsect)-min(lat_subsect),
max(lon_subsect)-min(lon_subsect)))
In [7]:
from utils.data_cube_utilities.dc_display_map import display_map
display_map(latitude = lat_subsect,longitude = lon_subsect)
Out[7]:
In [8]:
%%time
ds = dc.load(lat = lat_subsect,
lon = lon_subsect,
platform = platform,
product = product,
output_crs = "EPSG:32756",
measurements = ["red","blue","green","nir","quality"],
resolution = (-30,30))
In [9]:
ds
Out[9]:
In [10]:
from utils.data_cube_utilities.dc_rgb import rgb
rgb(ds.isel(time=6), x_coord='x', y_coord='y')
plt.show()
- The bathymetry function is located at the top of this notebook.
In [11]:
# Create Bathemtry Index column
ds["bathymetry"] = bathymetry_index(ds)
In [12]:
from utils.data_cube_utilities.dc_water_classifier import NDWI
# (green - nir) / (green + nir)
ds["ndwi"] = NDWI(ds, band_pair=1)
In [13]:
ds
Out[13]:
In [14]:
import os
from utils.data_cube_utilities.import_export import export_xarray_to_multiple_geotiffs
unmasked_dir = "geotiffs/landsat8/unmasked"
if not os.path.exists(unmasked_dir):
os.makedirs(unmasked_dir)
export_xarray_to_multiple_geotiffs(ds, unmasked_dir + "/unmasked.tif",
x_coord='x', y_coord='y')
In [15]:
# preview values
np.unique(ds["quality"])
Out[15]:
The threshold can be tuned if need be to better fit the RGB image above.
Unfortunately our existing WOFS algorithm is designed to work with Surface Reflectance (SR) and does not work with this data yet but with a few modifications it could be made to do so. We will approximate the WOFs mask withNDWIfor now.
In [16]:
# Tunable threshold for masking the land out
threshold = .05
water = (ds.ndwi>threshold).values
In [17]:
#preview one time slice to determine the effectiveness of the NDWI masking
rgb(ds.where(water).isel(time=6), x_coord='x', y_coord='y')
plt.show()
In [18]:
from utils.data_cube_utilities.dc_mosaic import ls8_oli_unpack_qa
clear_xarray = ls8_oli_unpack_qa(ds.quality, "clear")
In [19]:
full_mask = np.logical_and(clear_xarray, water)
ds = ds.where(full_mask)
In [20]:
plt.figure(figsize=[15,5])
#Visualize the distribution of the remaining data
sns.boxplot(ds['bathymetry'])
plt.show()
Interpretation: We can see that most of the values fall within a very short range. We can scale our plot's cmap limits to fit the specific quantile ranges for the bathymetry index so we can achieve better contrast from our plots.
In [21]:
#set the quantile range in either direction from the median value
def get_quantile_range(col, quantile_range = .25):
low = ds[col].quantile(.5 - quantile_range,["time","y","x"]).values
high = ds[col].quantile(.5 + quantile_range,["time","y","x"]).values
return low,high
In [22]:
#Custom function for a color mapping object
from matplotlib.colors import LinearSegmentedColormap
def custom_color_mapper(name = "custom", val_range = (1.96,1.96), colors = "RdGnBu"):
custom_cmap = LinearSegmentedColormap.from_list(name,colors=colors)
min, max = val_range
step = max/10.0
Z = [min,0],[0,max]
levels = np.arange(min,max+step,step)
cust_map = plt.contourf(Z, 100, cmap=custom_cmap)
plt.clf()
return cust_map.cmap
In [23]:
def mean_value_visual(ds, col, figsize = [15,15], cmap = "GnBu", low=None, high=None):
if low is None: low = np.min(ds[col]).values
if high is None: high = np.max(ds[col]).values
ds.reduce(np.nanmean,dim=["time"])[col].plot.imshow(figsize = figsize, cmap=cmap,
vmin=low, vmax=high)
In [24]:
mean_value_visual(ds, "bathymetry", cmap="GnBu")
If we clamp the range of the plot using different quantile ranges we can see relative differences in higher contrast.
In [25]:
# create range using the 10th and 90th quantile
low, high = get_quantile_range("bathymetry", .40)
custom = custom_color_mapper(val_range=(low,high),
colors=["darkred","red","orange","yellow","green",
"blue","darkblue","black"])
mean_value_visual(ds, "bathymetry", cmap=custom, low=low, high=high)
In [26]:
# create range using the 5th and 95th quantile
low, high = get_quantile_range("bathymetry", .45)
custom = custom_color_mapper(val_range=(low,high),
colors=["darkred","red","orange","yellow","green",
"blue","darkblue","black"])
mean_value_visual(ds, "bathymetry", cmap = custom, low=low, high = high)
In [27]:
# create range using the 2nd and 98th quantile
low, high = get_quantile_range("bathymetry", .48)
custom = custom_color_mapper(val_range=(low,high),
colors=["darkred","red","orange","yellow","green",
"blue","darkblue","black"])
mean_value_visual(ds, "bathymetry", cmap=custom, low=low, high=high)
In [28]:
# create range using the 1st and 99th quantile
low, high = get_quantile_range("bathymetry", .49)
custom = custom_color_mapper(val_range=(low,high),
colors=["darkred","red","orange","yellow","green",
"blue","darkblue","black"])
mean_value_visual(ds, "bathymetry", cmap=custom, low=low, high=high)
In [29]:
masked_dir = "geotiffs/landsat8/masked"
if not os.path.exists(masked_dir):
os.makedirs(masked_dir)
export_xarray_to_geotiff(ds, masked_dir + "/masked.tif",
x_coord='x', y_coord='y')