今天给大家推荐一篇高大上的文章:基于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
在当今这个数据驱动的时代,几乎每一个业务决策都离不开对数据的深入分析。而其中,指标波动归因分析更是至关重要的一环。无论是 ...
2025-02-18当数据开始说谎:那些年我们交过的学费 你有没有经历过这样的场景?熬了三个通宵做的数据分析报告,在会议上被老板一句"这数据靠 ...
2025-02-17数据分析作为一门跨学科领域,融合了统计学、编程、业务理解和可视化技术。无论是初学者还是有一定经验的从业者,系统化的学习路 ...
2025-02-17挖掘用户价值本质是让企业从‘赚今天的钱’升级为‘赚未来的钱’,同时让用户从‘被推销’变为‘被满足’。询问deepseek关于挖 ...
2025-02-17近来deepseek爆火,看看deepseek能否帮我们快速实现数据看板实时更新。 可以看出这对不知道怎么动手的小白来说是相当友好的, ...
2025-02-14一秒精通 Deepseek,不用找教程,不用买资料,更不用报一堆垃圾课程,所有这么去做的,都是舍近求远,因为你忽略了 deepseek 的 ...
2025-02-12自学 Python 的关键在于高效规划 + 实践驱动。以下是一份适合零基础快速入门的自学路径,结合资源推荐和实用技巧: 一、快速入 ...
2025-02-12“我们的利润率上升了,但销售额却没变,这是为什么?” “某个业务的市场份额在下滑,到底是什么原因?” “公司整体业绩 ...
2025-02-08活动介绍 为了助力大家在数据分析领域不断精进技能,我们特别举办本期打卡活动。在这里,你可以充分利用碎片化时间在线学习,让 ...
2025-02-071、闺女,醒醒,媒人把相亲的带来了。 我。。。。。。。 2、前年春节相亲相了40个, 去年春节相亲50个, 祖宗,今年你想相多少个 ...
2025-02-06在数据科学的广阔领域中,统计分析与数据挖掘占据了重要位置。尽管它们常常被视为有关联的领域,但两者在理论基础、目标、方法及 ...
2025-02-05在数据分析的世界里,“对比”是一种简单且有效的方法。这就像两个女孩子穿同一款式的衣服,效果不一样。 很多人都听过“货比三 ...
2025-02-05当我们只有非常少量的已标记数据,同时有大量未标记数据点时,可以使用半监督学习算法来处理。在sklearn中,基于图算法的半监督 ...
2025-02-05考虑一种棘手的情况:训练数据中大部分样本没有标签。此时,我们可以考虑使用半监督学习方法来处理。半监督学习能够利用这些额 ...
2025-02-04一、数学函数 1、取整 =INT(数字) 2、求余数 =MOD(除数,被除数) 3、四舍五入 =ROUND(数字,保留小数位数) 4、取绝对值 =AB ...
2025-02-03作者:CDA持证人 余治国 一般各平台出薪资报告,都会哀嚎遍野。举个例子,去年某招聘平台发布《中国女性职场现状调查报告》, ...
2025-02-02真正的数据分析大神是什么样的呢?有人认为他们能轻松驾驭各种分析工具,能够从海量数据中找到潜在关联,或者一眼识别报告中的数 ...
2025-02-01现今社会,“转行”似乎成无数职场人无法回避的话题。但行业就像座围城:外行人看光鲜,内行人看心酸。数据分析这个行业,近几年 ...
2025-01-31本人基本情况: 学校及专业:厦门大学经济学院应用统计 实习经历:快手数据分析、字节数据分析、百度数据分析 Offer情况:北京 ...
2025-01-3001专家简介 徐杨老师,CDA数据科学研究院教研副总监,主要负责CDA认证项目以及机器学习/人工智能类课程的研发与授课,负责过中 ...
2025-01-29