基于OpenCV实现海岸线变化检测

from glob import globimport numpy as np
import rasterioimport json, re, itertools, os
import matplotlib.pyplot as plt
import cv2 as cvfrom sklearn import preprocessingfrom sklearn.cluster import KMeans

N_OPTICS_BANDS = 7
with open("bands.json","r") as bandsJson:    bandsCharacteristics = json.load(bandsJson)

# bands.json
[{'id': '1', 'name': 'Coastal aerosol', 'span': '0.43-0.45', 'resolution': '30'}, {'id': '2', 'name': 'Blue', 'span': '0.45-0.51', 'resolution': '30'}, {'id': '3', 'name': 'Green', 'span': '0.53-0.59', 'resolution': '30'}, {'id': '4', 'name': 'Red', 'span': '0.64-0.67', 'resolution': '30'}, {'id': '5', 'name': 'NIR', 'span': '0.85-0.88', 'resolution': '30'}, {'id': '6', 'name': 'SWIR 1', 'span': '1.57-1.65', 'resolution': '30'}, {'id': '7', 'name': 'SWIR 2', 'span': '2.11-2.29', 'resolution': '30'}, {'id': '8', 'name': 'Panchromatic', 'span': '0.50-0.68', 'resolution': '15'}, {'id': '9', 'name': 'Cirrus', 'span': '1-36-1.38', 'resolution': '30'}, {'id': '10', 'name': 'TIRS 1', 'span': '10.6-11.9', 'resolution': '100'}, {'id': '11', 'name': 'TIRS 2', 'span': '11.50-12.51', 'resolution': '100'}]

class Acquisition:
def __init__(self, path, ext, nOpticsBands):self.nOpticsBands = nOpticsBandsself._getGeoTIFFs(path, ext)self.images = self._loadAcquisition()
def _getGeoTIFFs(self, path, ext):# It searches for GeoTIFF files within the folder.        print("Searching for '%s' files in %s" % (ext, path))self.fileList = glob(os.path.join(path,"*."+ext))self.opticsFileList = [    [list(filter(re.compile(r"band%s\."%a).search, self.fileList))[0] for a in range(1,self.nOpticsBands+1)]        print("Found %d 'tif' files" % len(self.opticsFileList))
def _loadAcquisition(self):# It finally reads and loads selected images into arrays.        print("Loading images")self.loads = [rasterio.open(bandPath) for bandPath in self.opticsFileList]        images = [load.read()[0] for load in self.loads]        print("Done")return images
def subsetImages(self, w1, w2, h1, h2, leftBound):# This function subsets images according the defined sizes.        print("Subsetting images (%s:%s, %s:%s)" % (w1, w2, h1, h2))        cols = (self.loads[0].bounds.left - leftBound)/30        registered = [np.insert(band,np.repeat(0,cols),0,axis=1) for band in self.images]        subset = [band[w1:w2,h1:h2] for band in registered]        print("Done")return subset

DATES = ["2014-02-01", "2019-07-25"]acquisitionsObjects = []
for date in DATES:    singleAcquisitionObject = Acquisition("Data/"+date, "tif", N_OPTICS_BANDS)    acquisitionsObjects.append( singleAcquisitionObject )

W1, W2 = 950, 2300H1, H2 = 4500, 5300
subAcquisitions = [acquisition.subsetImages(W1, W2, H1, H2, 552285.0) for acquisition in acquisitionsObjects].

axs = range(N_OPTICS_BANDS)fig, axs = plt.subplots(2, 4, figsize=(15,12))axs = list(itertools.chain.from_iterable(axs))
for b in range(N_OPTICS_BANDS):    id_ = bandsCharacteristics[b]["id"]    name_ = bandsCharacteristics[b]["name"]    span_ = bandsCharacteristics[b]["span"]    resolution_ = bandsCharacteristics[b]["resolution"]    title = "%s - %s\n%s (%s m)" % (id_, name_, span_, resolution_)    axs[b].imshow(preprocessing.StandardScaler().fit_transform(subAcquisitions[1][b]), cmap="Greys_r")    axs[b].set_title(title); axs[b].set_xticklabels([]); axs[b].set_yticklabels([])
plt.axis("off"); plt.tight_layout(w_pad=-10); plt.show()

BIAS = 1.5GAIN = [2.3,2.4,1.4]
r1 = (subAcquisitions[0][3] - subAcquisitions[0][3].min()) / (subAcquisitions[0][3].max()-subAcquisitions[0][3].min()) * GAIN[0] * BIASg1 = (subAcquisitions[0][2] - subAcquisitions[0][2].min()) / (subAcquisitions[0][2].max()-subAcquisitions[0][2].min()) * GAIN[1] * BIASb1 = (subAcquisitions[0][1] - subAcquisitions[0][1].min()) / (subAcquisitions[0][1].max()-subAcquisitions[0][1].min()) * GAIN[2] * BIAS
r2 = (subAcquisitions[1][3] - subAcquisitions[1][3].min()) / (subAcquisitions[1][3].max()-subAcquisitions[1][3].min()) * GAIN[0] * BIASg2 = (subAcquisitions[1][2] - subAcquisitions[1][2].min()) / (subAcquisitions[1][2].max()-subAcquisitions[1][2].min()) * GAIN[1] * BIASb2 = (subAcquisitions[1][1] - subAcquisitions[1][1].min()) / (subAcquisitions[1][1].max()-subAcquisitions[1][1].min()) * GAIN[2] * BIAS
rgbImage1, rgbImage2 = np.zeros((W2-W1,H2-H1,3)), np.zeros((W2-W1,H2-H1,3))rgbImage1[:,:,0], rgbImage2[:,:,0] = r1, r2rgbImage1[:,:,1], rgbImage2[:,:,1] = g1, g2rgbImage1[:,:,2], rgbImage2[:,:,2] = b1, b2
fig, (ax1,ax2) = plt.subplots(1,2,figsize=(16,12))ax1.imshow(rgbImage1); ax2.imshow(rgbImage2)ax1.set_title("RGB\n(Bands 4-3-2)\n2014-02-01"); ax2.set_title("RGB\n(Bands 4-3-2)\n2019-07-25")plt.show()

clusteredImages = [clusterLabels.reshape(subAcquisitions[0][0].shape).astype("uint8") for clusterLabels in clusters]blurredImages = [cv.GaussianBlur(clusteredImage, (15,15), 0) for clusteredImage in clusteredImages]
fig, (ax1, ax2) = plt.subplots(1,2,figsize=(16,13))
ax1.imshow(blurredImages[0])ax1.set_title("2014-02-01\nGaussian Blurred Image")ax2.imshow(blurredImages[1])ax2.set_title("2019-07-25\nGaussian Blurred Image")
plt.show()

rawEdges = [cv.Canny(blurredImage, 2, 5).astype("float").reshape(clusteredImages[0].shape) for blurredImage in blurredImages]
edges = []for edge in rawEdges:    edge[edge == 0] = np.nan    edges.append(edge)

plt.figure(figsize=(16,30))plt.imshow(rgbImage2)plt.imshow(edges[0], cmap = 'Set3_r')plt.imshow(edges[1], cmap = 'Set1')plt.title('CoastLine')plt.show()