forked from kilianv/CoastSat_WRL
You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
667 lines
27 KiB
Python
667 lines
27 KiB
Python
6 years ago
|
"""This module contains all the functions needed to preprocess the satellite images: creating a
|
||
|
|
cloud mask and pansharpening/downsampling the images.
|
||
|
|
|
||
|
|
Author: Kilian Vos, Water Research Laboratory, University of New South Wales
|
||
|
|
"""
|
||
|
|
|
||
|
|
# Initial settings
|
||
|
|
import os
|
||
|
|
import numpy as np
|
||
|
|
import matplotlib.pyplot as plt
|
||
|
|
from osgeo import gdal, ogr, osr
|
||
|
|
import skimage.transform as transform
|
||
|
|
import skimage.morphology as morphology
|
||
|
|
import sklearn.decomposition as decomposition
|
||
|
|
import skimage.exposure as exposure
|
||
|
|
from pylab import ginput
|
||
|
|
import pickle
|
||
|
|
import pdb
|
||
|
|
import SDS_tools
|
||
|
|
|
||
|
|
# Functions
|
||
|
|
|
||
|
|
def create_cloud_mask(im_qa, satname):
|
||
|
|
"""
|
||
|
|
Creates a cloud mask from the image containing the QA band information.
|
||
|
|
|
||
|
|
KV WRL 2018
|
||
|
|
|
||
|
|
Arguments:
|
||
|
|
-----------
|
||
|
|
im_qa: np.array
|
||
|
|
Image containing the QA band
|
||
|
|
satname: string
|
||
|
|
short name for the satellite (L8, L7, S2)
|
||
|
|
|
||
|
|
Returns:
|
||
|
|
-----------
|
||
|
|
cloud_mask : np.ndarray of booleans
|
||
|
|
A boolean array with True where the cloud are present
|
||
|
|
"""
|
||
|
|
|
||
|
|
# convert QA bits depending on the satellite mission
|
||
|
|
if satname == 'L8':
|
||
|
|
cloud_values = [2800, 2804, 2808, 2812, 6896, 6900, 6904, 6908]
|
||
|
|
elif satname == 'L7' or satname == 'L5' or satname == 'L4':
|
||
|
|
cloud_values = [752, 756, 760, 764]
|
||
|
|
elif satname == 'S2':
|
||
|
|
cloud_values = [1024, 2048] # 1024 = dense cloud, 2048 = cirrus clouds
|
||
|
|
|
||
|
|
# find which pixels have bits corresponding to cloud values
|
||
|
|
cloud_mask = np.isin(im_qa, cloud_values)
|
||
|
|
|
||
|
|
# remove isolated cloud pixels (there are some in the swash zone and they can cause problems)
|
||
|
|
if sum(sum(cloud_mask)) > 0 and sum(sum(~cloud_mask)) > 0:
|
||
|
|
morphology.remove_small_objects(cloud_mask, min_size=10, connectivity=1, in_place=True)
|
||
|
|
|
||
|
|
return cloud_mask
|
||
|
|
|
||
|
|
def hist_match(source, template):
|
||
|
|
"""
|
||
|
|
Adjust the pixel values of a grayscale image such that its histogram matches that of a
|
||
|
|
target image.
|
||
|
|
|
||
|
|
Arguments:
|
||
|
|
-----------
|
||
|
|
source: np.array
|
||
|
|
Image to transform; the histogram is computed over the flattened
|
||
|
|
array
|
||
|
|
template: np.array
|
||
|
|
Template image; can have different dimensions to source
|
||
|
|
Returns:
|
||
|
|
-----------
|
||
|
|
matched: np.array
|
||
|
|
The transformed output image
|
||
|
|
"""
|
||
|
|
|
||
|
|
oldshape = source.shape
|
||
|
|
source = source.ravel()
|
||
|
|
template = template.ravel()
|
||
|
|
|
||
|
|
# get the set of unique pixel values and their corresponding indices and
|
||
|
|
# counts
|
||
|
|
s_values, bin_idx, s_counts = np.unique(source, return_inverse=True,
|
||
|
|
return_counts=True)
|
||
|
|
t_values, t_counts = np.unique(template, return_counts=True)
|
||
|
|
|
||
|
|
# take the cumsum of the counts and normalize by the number of pixels to
|
||
|
|
# get the empirical cumulative distribution functions for the source and
|
||
|
|
# template images (maps pixel value --> quantile)
|
||
|
|
s_quantiles = np.cumsum(s_counts).astype(np.float64)
|
||
|
|
s_quantiles /= s_quantiles[-1]
|
||
|
|
t_quantiles = np.cumsum(t_counts).astype(np.float64)
|
||
|
|
t_quantiles /= t_quantiles[-1]
|
||
|
|
|
||
|
|
# interpolate linearly to find the pixel values in the template image
|
||
|
|
# that correspond most closely to the quantiles in the source image
|
||
|
|
interp_t_values = np.interp(s_quantiles, t_quantiles, t_values)
|
||
|
|
|
||
|
|
return interp_t_values[bin_idx].reshape(oldshape)
|
||
|
|
|
||
|
|
def pansharpen(im_ms, im_pan, cloud_mask):
|
||
|
|
"""
|
||
|
|
Pansharpens a multispectral image (3D), using the panchromatic band (2D) and a cloud mask.
|
||
|
|
A PCA is applied to the image, then the 1st PC is replaced with the panchromatic band.
|
||
|
|
|
||
|
|
KV WRL 2018
|
||
|
|
|
||
|
|
Arguments:
|
||
|
|
-----------
|
||
|
|
im_ms: np.array
|
||
|
|
Multispectral image to pansharpen (3D)
|
||
|
|
im_pan: np.array
|
||
|
|
Panchromatic band (2D)
|
||
|
|
cloud_mask: np.array
|
||
|
|
2D cloud mask with True where cloud pixels are
|
||
|
|
|
||
|
|
Returns:
|
||
|
|
-----------
|
||
|
|
im_ms_ps: np.ndarray
|
||
|
|
Pansharpened multisoectral image (3D)
|
||
|
|
"""
|
||
|
|
|
||
|
|
# reshape image into vector and apply cloud mask
|
||
|
|
vec = im_ms.reshape(im_ms.shape[0] * im_ms.shape[1], im_ms.shape[2])
|
||
|
|
vec_mask = cloud_mask.reshape(im_ms.shape[0] * im_ms.shape[1])
|
||
|
|
vec = vec[~vec_mask, :]
|
||
|
|
# apply PCA to RGB bands
|
||
|
|
pca = decomposition.PCA()
|
||
|
|
vec_pcs = pca.fit_transform(vec)
|
||
|
|
|
||
|
|
# replace 1st PC with pan band (after matching histograms)
|
||
|
|
vec_pan = im_pan.reshape(im_pan.shape[0] * im_pan.shape[1])
|
||
|
|
vec_pan = vec_pan[~vec_mask]
|
||
|
|
vec_pcs[:,0] = hist_match(vec_pan, vec_pcs[:,0])
|
||
|
|
vec_ms_ps = pca.inverse_transform(vec_pcs)
|
||
|
|
|
||
|
|
# reshape vector into image
|
||
|
|
vec_ms_ps_full = np.ones((len(vec_mask), im_ms.shape[2])) * np.nan
|
||
|
|
vec_ms_ps_full[~vec_mask,:] = vec_ms_ps
|
||
|
|
im_ms_ps = vec_ms_ps_full.reshape(im_ms.shape[0], im_ms.shape[1], im_ms.shape[2])
|
||
|
|
|
||
|
|
return im_ms_ps
|
||
|
|
|
||
|
|
|
||
|
|
def rescale_image_intensity(im, cloud_mask, prob_high):
|
||
|
|
"""
|
||
|
|
Rescales the intensity of an image (multispectral or single band) by applying
|
||
|
|
a cloud mask and clipping the prob_high upper percentile. This functions allows
|
||
|
|
to stretch the contrast of an image for visualisation purposes.
|
||
|
|
|
||
|
|
KV WRL 2018
|
||
|
|
|
||
|
|
Arguments:
|
||
|
|
-----------
|
||
|
|
im: np.array
|
||
|
|
Image to rescale, can be 3D (multispectral) or 2D (single band)
|
||
|
|
cloud_mask: np.array
|
||
|
|
2D cloud mask with True where cloud pixels are
|
||
|
|
prob_high: float
|
||
|
|
probability of exceedence used to calculate the upper percentile
|
||
|
|
|
||
|
|
Returns:
|
||
|
|
-----------
|
||
|
|
im_adj: np.array
|
||
|
|
The rescaled image
|
||
|
|
"""
|
||
|
|
|
||
|
|
# lower percentile is set to 0
|
||
|
|
prc_low = 0
|
||
|
|
|
||
|
|
# reshape the 2D cloud mask into a 1D vector
|
||
|
|
vec_mask = cloud_mask.reshape(im.shape[0] * im.shape[1])
|
||
|
|
|
||
|
|
# if image contains several bands, stretch the contrast for each band
|
||
|
|
if len(im.shape) > 2:
|
||
|
|
# reshape into a vector
|
||
|
|
vec = im.reshape(im.shape[0] * im.shape[1], im.shape[2])
|
||
|
|
# initiliase with NaN values
|
||
|
|
vec_adj = np.ones((len(vec_mask), im.shape[2])) * np.nan
|
||
|
|
# loop through the bands
|
||
|
|
for i in range(im.shape[2]):
|
||
|
|
# find the higher percentile (based on prob)
|
||
|
|
prc_high = np.percentile(vec[~vec_mask, i], prob_high)
|
||
|
|
# clip the image around the 2 percentiles and rescale the contrast
|
||
|
|
vec_rescaled = exposure.rescale_intensity(vec[~vec_mask, i],
|
||
|
|
in_range=(prc_low, prc_high))
|
||
|
|
vec_adj[~vec_mask,i] = vec_rescaled
|
||
|
|
# reshape into image
|
||
|
|
im_adj = vec_adj.reshape(im.shape[0], im.shape[1], im.shape[2])
|
||
|
|
|
||
|
|
# if image only has 1 bands (grayscale image)
|
||
|
|
else:
|
||
|
|
vec = im.reshape(im.shape[0] * im.shape[1])
|
||
|
|
vec_adj = np.ones(len(vec_mask)) * np.nan
|
||
|
|
prc_high = np.percentile(vec[~vec_mask], prob_high)
|
||
|
|
vec_rescaled = exposure.rescale_intensity(vec[~vec_mask], in_range=(prc_low, prc_high))
|
||
|
|
vec_adj[~vec_mask] = vec_rescaled
|
||
|
|
im_adj = vec_adj.reshape(im.shape[0], im.shape[1])
|
||
|
|
|
||
|
|
return im_adj
|
||
|
|
|
||
|
|
def preprocess_single(fn, satname):
|
||
|
|
"""
|
||
|
|
Creates a cloud mask using the QA band and performs pansharpening/down-sampling of the image.
|
||
|
|
|
||
|
|
KV WRL 2018
|
||
|
|
|
||
|
|
Arguments:
|
||
|
|
-----------
|
||
|
|
fn: str or list of str
|
||
|
|
filename of the .TIF file containing the image
|
||
|
|
for L7, L8 and S2 there is a filename for the bands at different resolutions
|
||
|
|
satname: str
|
||
|
|
name of the satellite mission (e.g., 'L5')
|
||
|
|
|
||
|
|
Returns:
|
||
|
|
-----------
|
||
|
|
im_ms: np.array
|
||
|
|
3D array containing the pansharpened/down-sampled bands (B,G,R,NIR,SWIR1)
|
||
|
|
georef: np.array
|
||
|
|
vector of 6 elements [Xtr, Xscale, Xshear, Ytr, Yshear, Yscale] defining the
|
||
|
|
coordinates of the top-left pixel of the image
|
||
|
|
cloud_mask: np.array
|
||
|
|
2D cloud mask with True where cloud pixels are
|
||
|
|
|
||
|
|
"""
|
||
|
|
|
||
|
|
#=============================================================================================#
|
||
|
|
# L5 images
|
||
|
|
#=============================================================================================#
|
||
|
|
if satname == 'L5':
|
||
|
|
|
||
|
|
# read all bands
|
||
|
|
data = gdal.Open(fn, gdal.GA_ReadOnly)
|
||
|
|
georef = np.array(data.GetGeoTransform())
|
||
|
|
bands = [data.GetRasterBand(k + 1).ReadAsArray() for k in range(data.RasterCount)]
|
||
|
|
im_ms = np.stack(bands, 2)
|
||
|
|
|
||
|
|
# down-sample to 15 m (half of the original pixel size)
|
||
|
|
nrows = im_ms.shape[0]*2
|
||
|
|
ncols = im_ms.shape[1]*2
|
||
|
|
|
||
|
|
# create cloud mask
|
||
|
|
im_qa = im_ms[:,:,5]
|
||
|
|
im_ms = im_ms[:,:,:-1]
|
||
|
|
cloud_mask = create_cloud_mask(im_qa, satname)
|
||
|
|
|
||
|
|
# resize the image using bilinear interpolation (order 1)
|
||
|
|
im_ms = transform.resize(im_ms,(nrows, ncols), order=1, preserve_range=True,
|
||
|
|
mode='constant')
|
||
|
|
# resize the image using nearest neighbour interpolation (order 0)
|
||
|
|
cloud_mask = transform.resize(cloud_mask, (nrows, ncols), order=0, preserve_range=True,
|
||
|
|
mode='constant').astype('bool_')
|
||
|
|
|
||
|
|
# adjust georeferencing vector to the new image size
|
||
|
|
# scale becomes 15m and the origin is adjusted to the center of new top left pixel
|
||
|
|
georef[1] = 15
|
||
|
|
georef[5] = -15
|
||
|
|
georef[0] = georef[0] + 7.5
|
||
|
|
georef[3] = georef[3] - 7.5
|
||
|
|
|
||
|
|
# check if -inf or nan values on any band and add to cloud mask
|
||
|
|
for k in range(im_ms.shape[2]):
|
||
|
|
im_inf = np.isin(im_ms[:,:,k], -np.inf)
|
||
|
|
im_nan = np.isnan(im_ms[:,:,k])
|
||
|
|
cloud_mask = np.logical_or(np.logical_or(cloud_mask, im_inf), im_nan)
|
||
|
|
|
||
|
|
# calculate cloud cover
|
||
|
|
cloud_cover = sum(sum(cloud_mask.astype(int)))/(cloud_mask.shape[0]*cloud_mask.shape[1])
|
||
|
|
|
||
|
|
#=============================================================================================#
|
||
|
|
# L7 images
|
||
|
|
#=============================================================================================#
|
||
|
|
elif satname == 'L7':
|
||
|
|
|
||
|
|
# read pan image
|
||
|
|
fn_pan = fn[0]
|
||
|
|
data = gdal.Open(fn_pan, gdal.GA_ReadOnly)
|
||
|
|
georef = np.array(data.GetGeoTransform())
|
||
|
|
bands = [data.GetRasterBand(k + 1).ReadAsArray() for k in range(data.RasterCount)]
|
||
|
|
im_pan = np.stack(bands, 2)[:,:,0]
|
||
|
|
|
||
|
|
# size of pan image
|
||
|
|
nrows = im_pan.shape[0]
|
||
|
|
ncols = im_pan.shape[1]
|
||
|
|
|
||
|
|
# read ms image
|
||
|
|
fn_ms = fn[1]
|
||
|
|
data = gdal.Open(fn_ms, gdal.GA_ReadOnly)
|
||
|
|
bands = [data.GetRasterBand(k + 1).ReadAsArray() for k in range(data.RasterCount)]
|
||
|
|
im_ms = np.stack(bands, 2)
|
||
|
|
|
||
|
|
# create cloud mask
|
||
|
|
im_qa = im_ms[:,:,5]
|
||
|
|
cloud_mask = create_cloud_mask(im_qa, satname)
|
||
|
|
|
||
|
|
# resize the image using bilinear interpolation (order 1)
|
||
|
|
im_ms = im_ms[:,:,:5]
|
||
|
|
im_ms = transform.resize(im_ms,(nrows, ncols), order=1, preserve_range=True,
|
||
|
|
mode='constant')
|
||
|
|
# resize the image using nearest neighbour interpolation (order 0)
|
||
|
|
cloud_mask = transform.resize(cloud_mask, (nrows, ncols), order=0, preserve_range=True,
|
||
|
|
mode='constant').astype('bool_')
|
||
|
|
# check if -inf or nan values on any band and eventually add those pixels to cloud mask
|
||
|
|
for k in range(im_ms.shape[2]+1):
|
||
|
|
if k == 5:
|
||
|
|
im_inf = np.isin(im_pan, -np.inf)
|
||
|
|
im_nan = np.isnan(im_pan)
|
||
|
|
else:
|
||
|
|
im_inf = np.isin(im_ms[:,:,k], -np.inf)
|
||
|
|
im_nan = np.isnan(im_ms[:,:,k])
|
||
|
|
cloud_mask = np.logical_or(np.logical_or(cloud_mask, im_inf), im_nan)
|
||
|
|
|
||
|
|
# calculate cloud cover
|
||
|
|
cloud_cover = sum(sum(cloud_mask.astype(int)))/(cloud_mask.shape[0]*cloud_mask.shape[1])
|
||
|
|
|
||
|
|
# pansharpen Green, Red, NIR (where there is overlapping with pan band in L7)
|
||
|
|
try:
|
||
|
|
im_ms_ps = pansharpen(im_ms[:,:,[1,2,3]], im_pan, cloud_mask)
|
||
|
|
except: # if pansharpening fails, keep downsampled bands (for long runs)
|
||
|
|
im_ms_ps = im_ms[:,:,[1,2,3]]
|
||
|
|
# add downsampled Blue and SWIR1 bands
|
||
|
|
im_ms_ps = np.append(im_ms[:,:,[0]], im_ms_ps, axis=2)
|
||
|
|
im_ms_ps = np.append(im_ms_ps, im_ms[:,:,[4]], axis=2)
|
||
|
|
|
||
|
|
im_ms = im_ms_ps.copy()
|
||
|
|
|
||
|
|
#=============================================================================================#
|
||
|
|
# L8 images
|
||
|
|
#=============================================================================================#
|
||
|
|
elif satname == 'L8':
|
||
|
|
|
||
|
|
# read pan image
|
||
|
|
fn_pan = fn[0]
|
||
|
|
data = gdal.Open(fn_pan, gdal.GA_ReadOnly)
|
||
|
|
georef = np.array(data.GetGeoTransform())
|
||
|
|
bands = [data.GetRasterBand(k + 1).ReadAsArray() for k in range(data.RasterCount)]
|
||
|
|
im_pan = np.stack(bands, 2)[:,:,0]
|
||
|
|
|
||
|
|
# size of pan image
|
||
|
|
nrows = im_pan.shape[0]
|
||
|
|
ncols = im_pan.shape[1]
|
||
|
|
|
||
|
|
# read ms image
|
||
|
|
fn_ms = fn[1]
|
||
|
|
data = gdal.Open(fn_ms, gdal.GA_ReadOnly)
|
||
|
|
bands = [data.GetRasterBand(k + 1).ReadAsArray() for k in range(data.RasterCount)]
|
||
|
|
im_ms = np.stack(bands, 2)
|
||
|
|
|
||
|
|
# create cloud mask
|
||
|
|
im_qa = im_ms[:,:,5]
|
||
|
|
cloud_mask = create_cloud_mask(im_qa, satname)
|
||
|
|
|
||
|
|
# resize the image using bilinear interpolation (order 1)
|
||
|
|
im_ms = im_ms[:,:,:5]
|
||
|
|
im_ms = transform.resize(im_ms,(nrows, ncols), order=1, preserve_range=True,
|
||
|
|
mode='constant')
|
||
|
|
# resize the image using nearest neighbour interpolation (order 0)
|
||
|
|
cloud_mask = transform.resize(cloud_mask, (nrows, ncols), order=0, preserve_range=True,
|
||
|
|
mode='constant').astype('bool_')
|
||
|
|
# check if -inf or nan values on any band and eventually add those pixels to cloud mask
|
||
|
|
for k in range(im_ms.shape[2]+1):
|
||
|
|
if k == 5:
|
||
|
|
im_inf = np.isin(im_pan, -np.inf)
|
||
|
|
im_nan = np.isnan(im_pan)
|
||
|
|
else:
|
||
|
|
im_inf = np.isin(im_ms[:,:,k], -np.inf)
|
||
|
|
im_nan = np.isnan(im_ms[:,:,k])
|
||
|
|
cloud_mask = np.logical_or(np.logical_or(cloud_mask, im_inf), im_nan)
|
||
|
|
|
||
|
|
# calculate cloud cover
|
||
|
|
cloud_cover = sum(sum(cloud_mask.astype(int)))/(cloud_mask.shape[0]*cloud_mask.shape[1])
|
||
|
|
|
||
|
|
# pansharpen Blue, Green, Red (where there is overlapping with pan band in L8)
|
||
|
|
try:
|
||
|
|
im_ms_ps = pansharpen(im_ms[:,:,[0,1,2]], im_pan, cloud_mask)
|
||
|
|
except: # if pansharpening fails, keep downsampled bands (for long runs)
|
||
|
|
im_ms_ps = im_ms[:,:,[0,1,2]]
|
||
|
|
# add downsampled NIR and SWIR1 bands
|
||
|
|
im_ms_ps = np.append(im_ms_ps, im_ms[:,:,[3,4]], axis=2)
|
||
|
|
|
||
|
|
im_ms = im_ms_ps.copy()
|
||
|
|
|
||
|
|
#=============================================================================================#
|
||
|
|
# S2 images
|
||
|
|
#=============================================================================================#
|
||
|
|
if satname == 'S2':
|
||
|
|
|
||
|
|
# read 10m bands (R,G,B,NIR)
|
||
|
|
fn10 = fn[0]
|
||
|
|
data = gdal.Open(fn10, gdal.GA_ReadOnly)
|
||
|
|
georef = np.array(data.GetGeoTransform())
|
||
|
|
bands = [data.GetRasterBand(k + 1).ReadAsArray() for k in range(data.RasterCount)]
|
||
|
|
im10 = np.stack(bands, 2)
|
||
|
|
im10 = im10/10000 # TOA scaled to 10000
|
||
|
|
|
||
|
|
# if image contains only zeros (can happen with S2), skip the image
|
||
|
|
if sum(sum(sum(im10))) < 1:
|
||
|
|
im_ms = []
|
||
|
|
georef = []
|
||
|
|
# skip the image by giving it a full cloud_mask
|
||
|
|
cloud_mask = np.ones((im10.shape[0],im10.shape[1])).astype('bool')
|
||
|
|
return im_ms, georef, cloud_mask
|
||
|
|
|
||
|
|
# size of 10m bands
|
||
|
|
nrows = im10.shape[0]
|
||
|
|
ncols = im10.shape[1]
|
||
|
|
|
||
|
|
# read 20m band (SWIR1)
|
||
|
|
fn20 = fn[1]
|
||
|
|
data = gdal.Open(fn20, gdal.GA_ReadOnly)
|
||
|
|
bands = [data.GetRasterBand(k + 1).ReadAsArray() for k in range(data.RasterCount)]
|
||
|
|
im20 = np.stack(bands, 2)
|
||
|
|
im20 = im20[:,:,0]
|
||
|
|
im20 = im20/10000 # TOA scaled to 10000
|
||
|
|
|
||
|
|
# resize the image using bilinear interpolation (order 1)
|
||
|
|
im_swir = transform.resize(im20, (nrows, ncols), order=1, preserve_range=True,
|
||
|
|
mode='constant')
|
||
|
|
im_swir = np.expand_dims(im_swir, axis=2)
|
||
|
|
|
||
|
|
# append down-sampled SWIR1 band to the other 10m bands
|
||
|
|
im_ms = np.append(im10, im_swir, axis=2)
|
||
|
|
|
||
|
|
# create cloud mask using 60m QA band (not as good as Landsat cloud cover)
|
||
|
|
fn60 = fn[2]
|
||
|
|
data = gdal.Open(fn60, gdal.GA_ReadOnly)
|
||
|
|
bands = [data.GetRasterBand(k + 1).ReadAsArray() for k in range(data.RasterCount)]
|
||
|
|
im60 = np.stack(bands, 2)
|
||
|
|
im_qa = im60[:,:,0]
|
||
|
|
cloud_mask = create_cloud_mask(im_qa, satname)
|
||
|
|
# resize the cloud mask using nearest neighbour interpolation (order 0)
|
||
|
|
cloud_mask = transform.resize(cloud_mask,(nrows, ncols), order=0, preserve_range=True,
|
||
|
|
mode='constant')
|
||
|
|
# check if -inf or nan values on any band and add to cloud mask
|
||
|
|
for k in range(im_ms.shape[2]):
|
||
|
|
im_inf = np.isin(im_ms[:,:,k], -np.inf)
|
||
|
|
im_nan = np.isnan(im_ms[:,:,k])
|
||
|
|
cloud_mask = np.logical_or(np.logical_or(cloud_mask, im_inf), im_nan)
|
||
|
|
|
||
|
|
# calculate cloud cover
|
||
|
|
cloud_cover = sum(sum(cloud_mask.astype(int)))/(cloud_mask.shape[0]*cloud_mask.shape[1])
|
||
|
|
|
||
|
|
return im_ms, georef, cloud_mask
|
||
|
|
|
||
|
|
|
||
|
|
def create_jpg(im_ms, cloud_mask, date, satname, filepath):
|
||
|
|
"""
|
||
|
|
Saves a .jpg file with the RGB image as well as the NIR and SWIR1 grayscale images.
|
||
|
|
|
||
|
|
KV WRL 2018
|
||
|
|
|
||
|
|
Arguments:
|
||
|
|
-----------
|
||
|
|
im_ms: np.array
|
||
|
|
3D array containing the pansharpened/down-sampled bands (B,G,R,NIR,SWIR1)
|
||
|
|
cloud_mask: np.array
|
||
|
|
2D cloud mask with True where cloud pixels are
|
||
|
|
date: str
|
||
|
|
String containing the date at which the image was acquired
|
||
|
|
satname: str
|
||
|
|
name of the satellite mission (e.g., 'L5')
|
||
|
|
|
||
|
|
Returns:
|
||
|
|
-----------
|
||
|
|
Saves a .jpg image corresponding to the preprocessed satellite image
|
||
|
|
|
||
|
|
"""
|
||
|
|
|
||
|
|
# rescale image intensity for display purposes
|
||
|
|
im_RGB = rescale_image_intensity(im_ms[:,:,[2,1,0]], cloud_mask, 99.9)
|
||
|
|
im_NIR = rescale_image_intensity(im_ms[:,:,3], cloud_mask, 99.9)
|
||
|
|
im_SWIR = rescale_image_intensity(im_ms[:,:,4], cloud_mask, 99.9)
|
||
|
|
|
||
|
|
# make figure
|
||
|
|
fig = plt.figure()
|
||
|
|
fig.set_size_inches([18,9])
|
||
|
|
fig.set_tight_layout(True)
|
||
|
|
# RGB
|
||
|
|
plt.subplot(131)
|
||
|
|
plt.axis('off')
|
||
|
|
plt.imshow(im_RGB)
|
||
|
|
plt.title(date + ' ' + satname, fontsize=16)
|
||
|
|
# NIR
|
||
|
|
plt.subplot(132)
|
||
|
|
plt.axis('off')
|
||
|
|
plt.imshow(im_NIR, cmap='seismic')
|
||
|
|
plt.title('Near Infrared', fontsize=16)
|
||
|
|
# SWIR
|
||
|
|
plt.subplot(133)
|
||
|
|
plt.axis('off')
|
||
|
|
plt.imshow(im_SWIR, cmap='seismic')
|
||
|
|
plt.title('Short-wave Infrared', fontsize=16)
|
||
|
|
# save figure
|
||
|
|
plt.rcParams['savefig.jpeg_quality'] = 100
|
||
|
|
fig.savefig(os.path.join(filepath,
|
||
|
|
date + '_' + satname + '.jpg'), dpi=150)
|
||
|
|
plt.close()
|
||
|
|
|
||
|
|
|
||
|
|
def preprocess_all_images(metadata, settings):
|
||
|
|
"""
|
||
|
|
Saves a .jpg image for all the file contained in metadata.
|
||
|
|
|
||
|
|
KV WRL 2018
|
||
|
|
|
||
|
|
Arguments:
|
||
|
|
-----------
|
||
|
|
sitename: str
|
||
|
|
name of the site (and corresponding folder)
|
||
|
|
metadata: dict
|
||
|
|
contains all the information about the satellite images that were downloaded
|
||
|
|
cloud_thresh: float
|
||
|
|
maximum fraction of cloud cover allowed in the images
|
||
|
|
|
||
|
|
Returns:
|
||
|
|
-----------
|
||
|
|
Generates .jpg files for all the satellite images avaialble
|
||
|
|
|
||
|
|
"""
|
||
|
|
|
||
|
|
sitename = settings['sitename']
|
||
|
|
cloud_thresh = settings['cloud_thresh']
|
||
|
|
|
||
|
|
# create subfolder to store the jpg files
|
||
|
|
filepath_jpg = os.path.join(os.getcwd(), 'data', sitename, 'jpg_files', 'preprocessed')
|
||
|
|
try:
|
||
|
|
os.makedirs(filepath_jpg)
|
||
|
|
except:
|
||
|
|
print('')
|
||
|
|
|
||
|
|
# loop through satellite list
|
||
|
|
for satname in metadata.keys():
|
||
|
|
# access the images
|
||
|
|
if satname == 'L5':
|
||
|
|
# access downloaded Landsat 5 images
|
||
|
|
filepath = os.path.join(os.getcwd(), 'data', sitename, satname, '30m')
|
||
|
|
filenames = os.listdir(filepath)
|
||
|
|
elif satname == 'L7':
|
||
|
|
# access downloaded Landsat 7 images
|
||
|
|
filepath_pan = os.path.join(os.getcwd(), 'data', sitename, 'L7', 'pan')
|
||
|
|
filepath_ms = os.path.join(os.getcwd(), 'data', sitename, 'L7', 'ms')
|
||
|
|
filenames_pan = os.listdir(filepath_pan)
|
||
|
|
filenames_ms = os.listdir(filepath_ms)
|
||
|
|
if (not len(filenames_pan) == len(filenames_ms)):
|
||
|
|
raise 'error: not the same amount of files for pan and ms'
|
||
|
|
filepath = [filepath_pan, filepath_ms]
|
||
|
|
filenames = filenames_pan
|
||
|
|
elif satname == 'L8':
|
||
|
|
# access downloaded Landsat 7 images
|
||
|
|
filepath_pan = os.path.join(os.getcwd(), 'data', sitename, 'L8', 'pan')
|
||
|
|
filepath_ms = os.path.join(os.getcwd(), 'data', sitename, 'L8', 'ms')
|
||
|
|
filenames_pan = os.listdir(filepath_pan)
|
||
|
|
filenames_ms = os.listdir(filepath_ms)
|
||
|
|
if (not len(filenames_pan) == len(filenames_ms)):
|
||
|
|
raise 'error: not the same amount of files for pan and ms'
|
||
|
|
filepath = [filepath_pan, filepath_ms]
|
||
|
|
filenames = filenames_pan
|
||
|
|
elif satname == 'S2':
|
||
|
|
# access downloaded Sentinel 2 images
|
||
|
|
filepath10 = os.path.join(os.getcwd(), 'data', sitename, satname, '10m')
|
||
|
|
filenames10 = os.listdir(filepath10)
|
||
|
|
filepath20 = os.path.join(os.getcwd(), 'data', sitename, satname, '20m')
|
||
|
|
filenames20 = os.listdir(filepath20)
|
||
|
|
filepath60 = os.path.join(os.getcwd(), 'data', sitename, satname, '60m')
|
||
|
|
filenames60 = os.listdir(filepath60)
|
||
|
|
if (not len(filenames10) == len(filenames20)) or (not len(filenames20) == len(filenames60)):
|
||
|
|
raise 'error: not the same amount of files for 10, 20 and 60 m'
|
||
|
|
filepath = [filepath10, filepath20, filepath60]
|
||
|
|
filenames = filenames10
|
||
|
|
|
||
|
|
# loop through images
|
||
|
|
for i in range(len(filenames)):
|
||
|
|
# image filename
|
||
|
|
fn = SDS_tools.get_filenames(filenames[i],filepath, satname)
|
||
|
|
# preprocess image (cloud mask + pansharpening/downsampling)
|
||
|
|
im_ms, georef, cloud_mask = preprocess_single(fn, satname)
|
||
|
|
# calculate cloud cover
|
||
|
|
cloud_cover = np.divide(sum(sum(cloud_mask.astype(int))),
|
||
|
|
(cloud_mask.shape[0]*cloud_mask.shape[1]))
|
||
|
|
# skip image if cloud cover is above threshold
|
||
|
|
if cloud_cover > cloud_thresh:
|
||
|
|
continue
|
||
|
|
# save .jpg with date and satellite in the title
|
||
|
|
date = filenames[i][:10]
|
||
|
|
create_jpg(im_ms, cloud_mask, date, satname, filepath_jpg)
|
||
|
|
|
||
|
|
def get_reference_sl(metadata, settings):
|
||
|
|
|
||
|
|
sitename = settings['sitename']
|
||
|
|
|
||
|
|
# check if reference shoreline already exists
|
||
|
|
filepath = os.path.join(os.getcwd(), 'data', sitename)
|
||
|
|
filename = sitename + '_ref_sl.pkl'
|
||
|
|
if filename in os.listdir(filepath):
|
||
|
|
print('Reference shoreline already exists and was loaded')
|
||
|
|
with open(os.path.join(filepath, sitename + '_ref_sl.pkl'), 'rb') as f:
|
||
|
|
refsl = pickle.load(f)
|
||
|
|
return refsl
|
||
|
|
|
||
|
|
else:
|
||
|
|
satname = 'S2'
|
||
|
|
# access downloaded Sentinel 2 images
|
||
|
|
filepath10 = os.path.join(os.getcwd(), 'data', sitename, satname, '10m')
|
||
|
|
filenames10 = os.listdir(filepath10)
|
||
|
|
filepath20 = os.path.join(os.getcwd(), 'data', sitename, satname, '20m')
|
||
|
|
filenames20 = os.listdir(filepath20)
|
||
|
|
filepath60 = os.path.join(os.getcwd(), 'data', sitename, satname, '60m')
|
||
|
|
filenames60 = os.listdir(filepath60)
|
||
|
|
if (not len(filenames10) == len(filenames20)) or (not len(filenames20) == len(filenames60)):
|
||
|
|
raise 'error: not the same amount of files for 10, 20 and 60 m'
|
||
|
|
for i in range(len(filenames10)):
|
||
|
|
# image filename
|
||
|
|
fn = [os.path.join(filepath10, filenames10[i]),
|
||
|
|
os.path.join(filepath20, filenames20[i]),
|
||
|
|
os.path.join(filepath60, filenames60[i])]
|
||
|
|
# preprocess image (cloud mask + pansharpening/downsampling)
|
||
|
|
im_ms, georef, cloud_mask = preprocess_single(fn, satname)
|
||
|
|
# calculate cloud cover
|
||
|
|
cloud_cover = np.divide(sum(sum(cloud_mask.astype(int))),
|
||
|
|
(cloud_mask.shape[0]*cloud_mask.shape[1]))
|
||
|
|
# skip image if cloud cover is above threshold
|
||
|
|
if cloud_cover > settings['cloud_thresh']:
|
||
|
|
continue
|
||
|
|
# rescale image intensity for display purposes
|
||
|
|
im_RGB = rescale_image_intensity(im_ms[:,:,[2,1,0]], cloud_mask, 99.9)
|
||
|
|
# make figure
|
||
|
|
fig = plt.figure()
|
||
|
|
fig.set_size_inches([18,9])
|
||
|
|
fig.set_tight_layout(True)
|
||
|
|
# RGB
|
||
|
|
plt.axis('off')
|
||
|
|
plt.imshow(im_RGB)
|
||
|
|
plt.title('click <skip> if image is not clear enough to digitize the shoreline.\n' +
|
||
|
|
'Otherwise click on <keep> and start digitizing the shoreline.\n' +
|
||
|
|
'When finished digitizing the shoreline click on the scroll wheel ' +
|
||
|
|
'(middle click).', fontsize=14)
|
||
|
|
plt.text(0, 0.9, 'keep', size=16, ha="left", va="top",
|
||
|
|
transform=plt.gca().transAxes,
|
||
|
|
bbox=dict(boxstyle="square", ec='k',fc='w'))
|
||
|
|
plt.text(1, 0.9, 'skip', size=16, ha="right", va="top",
|
||
|
|
transform=plt.gca().transAxes,
|
||
|
|
bbox=dict(boxstyle="square", ec='k',fc='w'))
|
||
|
|
mng = plt.get_current_fig_manager()
|
||
|
|
mng.window.showMaximized()
|
||
|
|
# let user click on the image once
|
||
|
|
pt_keep = ginput(n=1, timeout=100, show_clicks=True)
|
||
|
|
pt_keep = np.array(pt_keep)
|
||
|
|
# if clicks next to <skip>, show another image
|
||
|
|
if pt_keep[0][0] > im_ms.shape[1]/2:
|
||
|
|
plt.close()
|
||
|
|
continue
|
||
|
|
else:
|
||
|
|
# let user click on the shoreline
|
||
|
|
pts = ginput(n=5000, timeout=100000, show_clicks=True)
|
||
|
|
pts_pix = np.array(pts)
|
||
|
|
plt.close()
|
||
|
|
# convert image coordinates to world coordinates
|
||
|
|
pts_world = SDS_tools.convert_pix2world(pts_pix[:,[1,0]], georef)
|
||
|
|
image_epsg = metadata[satname]['epsg'][i]
|
||
|
|
pts_coords = SDS_tools.convert_epsg(pts_world, image_epsg, settings['output_epsg'])
|
||
|
|
with open(os.path.join(filepath, sitename + '_ref_sl.pkl'), 'wb') as f:
|
||
|
|
pickle.dump(pts_coords, f)
|
||
|
|
print('Reference shoreline has been saved')
|
||
|
|
break
|
||
|
|
|
||
|
|
return pts_coords
|