今天給大家推薦一篇高大上的文章：基于OpenCV實現海岸線變化檢測。OpenCV大家都知道，一款開源的計算機視覺庫，平常大家看到的都是OpenCV人臉識別，圖像處理之類的，今天跟小編一起來看他如何實現海岸線變化檢測的吧。

文章來源： 小白學視覺

作者：努比

介紹

海岸是一個動態系統，其中因侵蝕現象導致的海岸線后退、或是由眾多因素如氣象，地質，生物和人類活動所導致線前進的是常見現象。

在海洋磨損作用大于沉積物的情況下，有明顯的海岸侵蝕，我們稱之為地球表面的崩解和破壞。

資料來源：弗林德斯大學（CC0）

本文的目標

在本文中，我們將對Landsat 8平臺上的OLI（陸地成像儀）傳感器獲取的衛星圖像使用Canny Edge Detection算法。

通過這種方法，我們將能夠可視化的估計特定歐洲地區遭受強腐蝕作用的海岸線隨時間的推移：霍德內斯海岸。

一下是處理流程：

處理流程

在開始之前讓我們先介紹一下OLI數據...

0.關于Landsat OLI數據的簡要介紹

Landsat 8是一個軌道平臺，安裝在稱為OLI（陸地成像儀）的11波段多光譜傳感器上。

具體來說，在本文中，我們將僅使用分辨率為30米（即前7個）的波段。

美國地質調查局陸地衛星8號

該數據可以免費下載，注冊后，獲得USGS：https://earthexplorer.usgs.gov/。

而且，通常我摸并不使用入射太陽光作為原始數據，而是使用反射率，即從地球表面反射的太陽光量[0-1]。

1.包導入

在各種常見的包，我們將使用rasterio處理圖像，利用OpenCV中的Canny 算法和Scikit-Learn分割圖像。

from glob import glob
import numpy as np

import rasterio
import json, re, itertools, os

import matplotlib.pyplot as plt

import cv2 as cv
from sklearn import preprocessing
from sklearn.cluster import KMeans

2.數據導入

讓我們定義一個變量，該變量告訴我們要保留的波段數以及在JSON中輸入的輔助數據：

N_OPTICS_BANDS = 7

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

這個Json是Landsat OLI成像儀的信息集合。類似于一種說明手冊：

# 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'}]

bands.json文件包含有關我們將要使用的頻段的所有有用信息。

注意，我們將僅使用分辨率為30 m的頻段，因此僅使用前7個頻段。如果您愿意使用較低的分辨率（100m），則也可以嵌入TIRS 1和TIRS 2頻段。

正如上面幾行已經提到的那樣，我們將使用從Landsat-8 OLI上獲取兩組不同的數據：

? 2014/02/01

? 2019/07/25

為了簡化兩次采集的所需操作，我們將定義一個Acquisition（）類，其中將封裝所有必要的函數。

在執行代碼期間，我們能夠執行一些基礎支持性的功能，例如：

? 在指定路徑中搜索GeoTIFF；

? 加載采購；

? 購置登記 （調整）；

? 收購子集

class Acquisition:
    
    def __init__(self, path, ext, nOpticsBands):
        self.nOpticsBands = nOpticsBands
        self._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 )

運行結果如下：

Searching for 'tif' files in Data/2014-02-01
Found 7 'tif' files
Loading images
Done
Searching for 'tif' files in Data/2019-07-25
Found 7 'tif' files
Loading images
Done

現在我們已加載了14張OLI圖像（在7個波段中各采集2個）。

2.1 子集多光譜立方體

在這個階段中，先對兩個多光譜立方體進行“對齊”（或正式注冊），再切出不感興趣的部分。

我們可以使用ImageImages（）函數“剪切”不需要的數據。

因此，我們定義AOI（感興趣的區域），并使用Acquisition（）類中的subsetImages（）函數進行設置：

W1, W2 = 950, 2300
H1, H2 = 4500, 5300

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

完成！

3.數據探索

3.1可視化多光譜立方體

讓我們嘗試查看2019/07/25收購的所有范圍。出于純粹的美學原因，在繪制圖像之前，讓我們使用

StandardScaler（）對圖像進行標準化。

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()

以下是運行結果。

這些圖中，有些波段比其他波段更亮。這很正常。

3.2可視化復合RGB中的多光譜立方體

現在，讓我們嘗試可視化使用波段4（紅色），3（綠色）和2（藍色）獲得的RGB復合圖像中的兩次采集。

定義BIAS和GAIN 僅是為了獲得更好的效果。

BIAS = 1.5
GAIN = [2.3,2.4,1.4]

r1 = (subAcquisitions[0][3] - subAcquisitions[0][3].min()) / (subAcquisitions[0][3].max()-subAcquisitions[0][3].min()) * GAIN[0] * BIAS
g1 = (subAcquisitions[0][2] - subAcquisitions[0][2].min()) / (subAcquisitions[0][2].max()-subAcquisitions[0][2].min()) * GAIN[1] * BIAS
b1 = (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] * BIAS
g2 = (subAcquisitions[1][2] - subAcquisitions[1][2].min()) / (subAcquisitions[1][2].max()-subAcquisitions[1][2].min()) * GAIN[1] * BIAS
b2 = (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, r2
rgbImage1[:,:,1], rgbImage2[:,:,1] = g1, g2
rgbImage1[:,:,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()

結果如下圖所示！有趣的是，這兩次獲取的反射率完全的不同。

好的，繼續進行海岸線檢測。

4.自動化海岸線檢測

在本段中，我們將使用Canny的算法執行邊緣檢測。

在進行實際檢測之前，有必要準備數據集，嘗試通過聚類算法對數據集進行分割以區分海洋和陸地。

4.1數據準備

在此階段，我們將重塑兩個多光譜立方體以進行聚類操作。

4.2用K均值進行圖像分割

我們通過k均值對這兩次采集進行細分（使用自己喜歡的模型即可）。

4.3細分結果

這是確定的代表新興土地和水體的兩個集群。

4.4Canny邊緣檢測算法

Canny的傳統鍵技術分為以下幾個階段：

1. 高斯濾波器通過卷積降低噪聲；

2. 四個方向（水平，垂直和2個傾斜）的圖像梯度計算；

3. 梯度局部最大值的提取;

4. 帶有滯后的閾值，用于邊緣提取。

讓我們開始，將聚類結果轉換為圖像，然后通過具有15x15內核的高斯濾波器降低噪聲：

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()

在圖像稍微模糊之后，我們可以使用OpenCV  Canny（）模塊：

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)

在單行代碼中，我們獲得了梯度，提取了局部最大值，然后對每次采集都應用了帶有滯后的閾值。

注意：我們可以使用不同參數Canny（）來探索處理結果。

4.5結果

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()

以下是一些詳細信息：

5結論

從結果中可以看到，Canny的算法在其原始管道中運行良好，但其性能通常取決于所涉及的數據。

實際上，所使用的聚類算法使我們能夠對多光譜立方體進行細分。并行使用多個聚類模型可以總體上改善結果。