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

文章來源： 小白學(xué)視覺

作者：努比

介紹

海岸是一個動態(tài)系統(tǒng)，其中因侵蝕現(xiàn)象導(dǎo)致的海岸線后退、或是由眾多因素如氣象，地質(zhì)，生物和人類活動所導(dǎo)致線前進(jìn)的是常見現(xiàn)象。

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

資料來源：弗林德斯大學(xué)（CC0）

本文的目標(biāo)

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

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

一下是處理流程：

處理流程

在開始之前讓我們先介紹一下OLI數(shù)據(jù)...

0.關(guān)于Landsat OLI數(shù)據(jù)的簡要介紹

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

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

美國地質(zhì)調(diào)查局陸地衛(wèi)星8號

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

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

1.包導(dǎo)入

在各種常見的包，我們將使用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.數(shù)據(jù)導(dǎo)入

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

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文件包含有關(guān)我們將要使用的頻段的所有有用信息。

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

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

? 2014/02/01

? 2019/07/25

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

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

? 在指定路徑中搜索GeoTIFF；

? 加載采購；

? 購置登記 （調(diào)整）；

? 收購子集

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

好的，讓我們現(xiàn)在開始啟動整個代碼：

DATES = ["2014-02-01", "2019-07-25"]
acquisitionsObjects = []

for date in DATES:
    singleAcquisitionObject = Acquisition("Data/"+date, "tif", N_OPTICS_BANDS)
    acquisitionsObjects.append( singleAcquisitionObject )

運行結(jié)果如下：

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

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

2.1 子集多光譜立方體

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

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

因此，我們定義AOI（感興趣的區(qū)域），并使用Acquisition（）類中的subsetImages（）函數(shù)進(jìn)行設(shè)置：

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

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

完成！

3.數(shù)據(jù)探索

3.1可視化多光譜立方體

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

StandardScaler（）對圖像進(jìn)行標(biāo)準(zhǔn)化。

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

以下是運行結(jié)果。

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

3.2可視化復(fù)合RGB中的多光譜立方體

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

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

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

好的，繼續(xù)進(jìn)行海岸線檢測。

4.自動化海岸線檢測

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

在進(jìn)行實際檢測之前，有必要準(zhǔn)備數(shù)據(jù)集，嘗試通過聚類算法對數(shù)據(jù)集進(jìn)行分割以區(qū)分海洋和陸地。

4.1數(shù)據(jù)準(zhǔn)備

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

4.2用K均值進(jìn)行圖像分割

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

4.3細(xì)分結(jié)果

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

4.4Canny邊緣檢測算法

Canny的傳統(tǒng)鍵技術(shù)分為以下幾個階段：

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

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

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

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

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

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)

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

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

4.5結(jié)果

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

以下是一些詳細(xì)信息：

5結(jié)論

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

實際上，所使用的聚類算法使我們能夠?qū)Χ喙庾V立方體進(jìn)行細(xì)分。并行使用多個聚類模型可以總體上改善結(jié)果。