
今天给大家推荐一篇高大上的文章:基于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的算法在其原始管道中运行良好,但其性能通常取决于所涉及的数据。
实际上,所使用的聚类算法使我们能够对多光谱立方体进行细分。并行使用多个聚类模型可以总体上改善结果。
数据分析咨询请扫描二维码
若不方便扫码,搜微信号:CDAshujufenxi
LSTM 模型输入长度选择技巧:提升序列建模效能的关键 在循环神经网络(RNN)家族中,长短期记忆网络(LSTM)凭借其解决长序列 ...
2025-07-11CDA 数据分析师报考条件详解与准备指南 在数据驱动决策的时代浪潮下,CDA 数据分析师认证愈发受到瞩目,成为众多有志投身数 ...
2025-07-11数据透视表中两列相乘合计的实用指南 在数据分析的日常工作中,数据透视表凭借其强大的数据汇总和分析功能,成为了 Excel 用户 ...
2025-07-11尊敬的考生: 您好! 我们诚挚通知您,CDA Level I和 Level II考试大纲将于 2025年7月25日 实施重大更新。 此次更新旨在确保认 ...
2025-07-10BI 大数据分析师:连接数据与业务的价值转化者 在大数据与商业智能(Business Intelligence,简称 BI)深度融合的时代,BI ...
2025-07-10SQL 在预测分析中的应用:从数据查询到趋势预判 在数据驱动决策的时代,预测分析作为挖掘数据潜在价值的核心手段,正被广泛 ...
2025-07-10数据查询结束后:分析师的收尾工作与价值深化 在数据分析的全流程中,“query end”(查询结束)并非工作的终点,而是将数 ...
2025-07-10CDA 数据分析师考试:从报考到取证的全攻略 在数字经济蓬勃发展的今天,数据分析师已成为各行业争抢的核心人才,而 CDA(Certi ...
2025-07-09【CDA干货】单样本趋势性检验:捕捉数据背后的时间轨迹 在数据分析的版图中,单样本趋势性检验如同一位耐心的侦探,专注于从单 ...
2025-07-09year_month数据类型:时间维度的精准切片 在数据的世界里,时间是最不可或缺的维度之一,而year_month数据类型就像一把精准 ...
2025-07-09CDA 备考干货:Python 在数据分析中的核心应用与实战技巧 在 CDA 数据分析师认证考试中,Python 作为数据处理与分析的核心 ...
2025-07-08SPSS 中的 Mann-Kendall 检验:数据趋势与突变分析的有力工具 在数据分析的广袤领域中,准确捕捉数据的趋势变化以及识别 ...
2025-07-08备战 CDA 数据分析师考试:需要多久?如何规划? CDA(Certified Data Analyst)数据分析师认证作为国内权威的数据分析能力认证 ...
2025-07-08LSTM 输出不确定的成因、影响与应对策略 长短期记忆网络(LSTM)作为循环神经网络(RNN)的一种变体,凭借独特的门控机制,在 ...
2025-07-07统计学方法在市场调研数据中的深度应用 市场调研是企业洞察市场动态、了解消费者需求的重要途径,而统计学方法则是市场调研数 ...
2025-07-07CDA数据分析师证书考试全攻略 在数字化浪潮席卷全球的当下,数据已成为企业决策、行业发展的核心驱动力,数据分析师也因此成为 ...
2025-07-07剖析 CDA 数据分析师考试题型:解锁高效备考与答题策略 CDA(Certified Data Analyst)数据分析师考试作为衡量数据专业能力的 ...
2025-07-04SQL Server 字符串截取转日期:解锁数据处理的关键技能 在数据处理与分析工作中,数据格式的规范性是保证后续分析准确性的基础 ...
2025-07-04CDA 数据分析师视角:从数据迷雾中探寻商业真相 在数字化浪潮席卷全球的今天,数据已成为企业决策的核心驱动力,CDA(Certifie ...
2025-07-04CDA 数据分析师:开启数据职业发展新征程 在数据成为核心生产要素的今天,数据分析师的职业价值愈发凸显。CDA(Certified D ...
2025-07-03