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

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

作者：努比

介紹

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

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

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

本文的目標(biāo)

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

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

一下是處理流程：

處理流程

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

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

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

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

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

該數(shù)據(jù)可以免費(fèi)下載，注冊(cè)后，獲得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)入

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

N_OPTICS_BANDS = 7

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

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

# 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個(gè)頻段。如果您愿意使用較低的分辨率（100m），則也可以嵌入TIRS 1和TIRS 2頻段。

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

? 2014/02/01

? 2019/07/25

為了簡(jiǎn)化兩次采集的所需操作，我們將定義一個(gè)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)在開始啟動(dòng)整個(gè)代碼：

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

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

運(yùn)行結(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個(gè)波段中各采集2個(gè)）。

2.1 子集多光譜立方體

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

我們可以使用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（）對(duì)圖像進(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()

以下是運(yùn)行結(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)行海岸線檢測(cè)。

4.自動(dòng)化海岸線檢測(cè)

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

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

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

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

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

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

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

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

4.4Canny邊緣檢測(cè)算法

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

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

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

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)

在單行代碼中，我們獲得了梯度，提取了局部最大值，然后對(duì)每次采集都應(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的算法在其原始管道中運(yùn)行良好，但其性能通常取決于所涉及的數(shù)據(jù)。

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