satellite coastline buffer.py

import cv2
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt


def buffer_coastline(img, threshold=135):
    arr = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)
    mask = np.ndarray(arr.shape + (3,)).astype(np.uint8)
    below = arr < threshold
    mask[below] = [255,255,255]
    mask = cv2.morphologyEx(mask[:,:,0], cv2.MORPH_OPEN, kernel=np.ones([3, 3]), iterations=3)
    edm1 = cv2.distanceTransform(mask, cv2.cv.CV_DIST_L2, 5)
    s = (edm1 < 10)
    mask[s] = 255
    mask[~s] = 0
    edm2 = cv2.distanceTransform(mask, cv2.cv.CV_DIST_L2, 5)
    s = (edm2 < 20) & (edm1 < 10)
    mask[s] = 255
    mask[~s] = 0
    return mask, edm1, edm2


def render_array(arr1, arr2=None, width=7, height=7, cmap=plt.cm.jet):
    if arr2 is None:
        fig, ax = plt.subplots(1, figsize=(width,height), facecolor='white')
        ax.axis('off')
        plt.imshow(arr1, cmap=cmap)
        plt.show()
    else: 
        fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2)
        fig.set_size_inches(width, height)
        ax1.imshow(arr2)
        ax1.axis('off')
        ax2.imshow(arr1)
        ax2.axis('off')
                          
            
img = clipped[3][:,:,0:3].copy()
mask, edm1, edm2 = buffer_coastline(img.copy())
render_array(edm1, edm2, width=10, height=10)

img = clipped[3].copy()
img[:,:,3] = mask
render_array(img, mask, width=10, height=10)

masks = []
for image_arr in clipped:
    mask, edm1, edm2 = buffer_coastline(image_arr, threshold=125)
    masks.append(mask)
    
union_mask = reduce( lambda s,e: s & e, masks )
union_mask = cv2.dilate(union_mask, kernel=np.ones([3, 3]), iterations=3)
binary_mask = union_mask.copy() 
binary_mask[(union_mask > 0)] = 1

img = clipped[3].copy()
img[:,:,3] = union_mask

render_array(img, width=10, height=10)