今天给大家推荐一篇高大上的文章:基于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
你是否被统计学复杂的理论和晦涩的公式劝退过?别担心,“山有木兮:统计学极简入门(Python)” 将为你一一化解这些难题。课程 ...
2025-03-31在电商、零售、甚至内容付费业务中,你真的了解你的客户吗? 有些客户下了一两次单就消失了,有些人每个月都回购,有些人曾经是 ...
2025-03-31在数字化浪潮中,数据驱动决策已成为企业发展的核心竞争力,数据分析人才的需求持续飙升。世界经济论坛发布的《未来就业报告》, ...
2025-03-28你有没有遇到过这样的情况?流量进来了,转化率却不高,辛辛苦苦拉来的用户,最后大部分都悄无声息地离开了,这时候漏斗分析就非 ...
2025-03-27TensorFlow Datasets(TFDS)是一个用于下载、管理和预处理机器学习数据集的库。它提供了易于使用的API,允许用户从现有集合中 ...
2025-03-26"不谋全局者,不足谋一域。"在数据驱动的商业时代,战略级数据分析能力已成为职场核心竞争力。《CDA二级教材:商业策略数据分析 ...
2025-03-26当你在某宝刷到【猜你喜欢】时,当抖音精准推来你的梦中情猫时,当美团外卖弹窗刚好是你想吃的火锅店…… 恭喜你,你正在被用户 ...
2025-03-26当面试官问起随机森林时,他到底在考察什么? ""请解释随机森林的原理""——这是数据分析岗位面试中的经典问题。但你可能不知道 ...
2025-03-25在数字化浪潮席卷的当下,数据俨然成为企业的命脉,贯穿于业务运作的各个环节。从线上到线下,从平台的交易数据,到门店的运营 ...
2025-03-25在互联网和移动应用领域,DAU(日活跃用户数)是一个耳熟能详的指标。无论是产品经理、运营,还是数据分析师,DAU都是衡量产品 ...
2025-03-24ABtest做的好,产品优化效果差不了!可见ABtest在评估优化策略的效果方面地位还是很高的,那么如何在业务中应用ABtest? 结合企业 ...
2025-03-21在企业数据分析中,指标体系是至关重要的工具。不仅帮助企业统一数据标准、提升数据质量,还能为业务决策提供有力支持。本文将围 ...
2025-03-20解锁数据分析师高薪密码,CDA 脱产就业班助你逆袭! 在数字化浪潮中,数据驱动决策已成为企业发展的核心竞争力,数据分析人才的 ...
2025-03-19在 MySQL 数据库中,查询一张表但是不包含某个字段可以通过以下两种方法实现:使用 SELECT 子句以明确指定想要的字段,或者使 ...
2025-03-17在当今数字化时代,数据成为企业发展的关键驱动力,而用户画像作为数据分析的重要成果,改变了企业理解用户、开展业务的方式。无 ...
2025-03-172025年是智能体(AI Agent)的元年,大模型和智能体的发展比较迅猛。感觉年初的deepseek刚火没多久,这几天Manus又成为媒体头条 ...
2025-03-14以下的文章内容来源于柯家媛老师的专栏,如果您想阅读专栏《小白必备的数据思维课》,点击下方链接 https://edu.cda.cn/goods/sh ...
2025-03-13以下的文章内容来源于刘静老师的专栏,如果您想阅读专栏《10大业务分析模型突破业务瓶颈》,点击下方链接 https://edu.cda.cn/go ...
2025-03-12以下的文章内容来源于柯家媛老师的专栏,如果您想阅读专栏《小白必备的数据思维课》,点击下方链接 https://edu.cda.cn/goods/sh ...
2025-03-11随着数字化转型的加速,企业积累了海量数据,如何从这些数据中挖掘有价值的信息,成为企业提升竞争力的关键。CDA认证考试体系应 ...
2025-03-10