python实现形态学建筑物指数MBI提取建筑物及数据获取

2023-12-13 09:28

本文主要是介绍python实现形态学建筑物指数MBI提取建筑物及数据获取,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

前言

    形态学建筑物指数MBI通过建立建筑物的隐式特征和形态学算子之间的关系进行建筑物的提取[1]。

原理

图片

上图源自[2]。

实验数据

简单找了一张小图片:

图片

test.jpg

代码

为了支持遥感图像,读写数据函数都是利用GDAL写的。

import numpy as np
import gdal#  读取tif数据集
def readTif(fileName, xoff = 0, yoff = 0, data_width = 0, data_height = 0):dataset = gdal.Open(fileName)if dataset == None:print(fileName + "文件无法打开")#  栅格矩阵的列数width = dataset.RasterXSize #  栅格矩阵的行数height = dataset.RasterYSize #  波段数bands = dataset.RasterCount #  获取数据if(data_width == 0 and data_height == 0):data_width = widthdata_height = heightdata = dataset.ReadAsArray(xoff, yoff, data_width, data_height)#  获取仿射矩阵信息geotrans = dataset.GetGeoTransform()#  获取投影信息proj = dataset.GetProjection()return width, height, bands, data, geotrans, proj#  保存tif文件函数
def writeTiff(im_data, im_geotrans, im_proj, path):if 'int8' in im_data.dtype.name:datatype = gdal.GDT_Byteelif 'int16' in im_data.dtype.name:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float32if len(im_data.shape) == 3:im_bands, im_height, im_width = im_data.shapeelif len(im_data.shape) == 2:im_data = np.array([im_data])im_bands, im_height, im_width = im_data.shape#创建文件driver = gdal.GetDriverByName("GTiff")dataset = driver.Create(path, int(im_width), int(im_height), int(im_bands), datatype)if(dataset!= None):dataset.SetGeoTransform(im_geotrans) #写入仿射变换参数dataset.SetProjection(im_proj) #写入投影for i in range(im_bands):dataset.GetRasterBand(i+1).WriteArray(im_data[i])del dataset

接下来就是就算MBI,代码注释很详细,也可以对着原理来看。

from skimage.morphology import square, white_tophat
from skimage.transform import rotate#  计算MBI
#  s_min: 结构元素大小最小值
#  s_max: 结构元素大小最大值
#  delta_s: 颗粒测定的间隔
def CalculationMBI(filePath, MBIPath, s_min, s_max, delta_s):#  读取图像的相关信息width, height, bands, image, geotrans, proj = readTif(filePath)#  多光谱带的最大值对应于具有高反射率的特征->取光谱带最大值作为后续计算数据gray = np.max(image, 0)#  为消除白帽边缘效应,进行边缘补零gray = np.pad(gray, ((s_min, s_min), (s_min, s_min)), 'constant', constant_values=(0, 0))#  形态学剖面集合MP_MBI_list = []#  差分形态学剖面DMP集合DMP_MBI_list = []#  计算形态学剖面for i in range(s_min, s_max + 1, 2 * delta_s):print("s = ", i)#  大小为i×i的单位矩阵SE_intermediate = square(i)#  只保留中间一行为1,其他设置为0SE_intermediate[ : int((i - 1) / 2), :] = 0SE_intermediate[int(((i - 1) / 2) + 1) : , :] = 0#  SE_intermediate表示结构元素,用于设定局部区域的形状和大小#  旋转0 45 90 135°for angle in range(0, 180, 45):SE_intermediate = rotate(SE_intermediate, angle, order = 0, preserve_range = True).astype('uint8')#  多角度形态学白帽重构MP_MBI = white_tophat(gray, selem = SE_intermediate)MP_MBI_list.append(MP_MBI)#  计算差分形态学剖面DMPfor j in range(4, len(MP_MBI_list), 1):#  差的绝对值DMP_MBI = np.absolute(MP_MBI_list[j] - MP_MBI_list[j - 4])DMP_MBI_list.append(DMP_MBI)#  计算MBIMBI = np.sum(DMP_MBI_list, axis = 0) / (4 * (((s_max - s_min) / delta_s) + 1))#  去除多余边缘结果MBI = MBI[s_min : MBI.shape[0] - s_min, s_min : MBI.shape[1] - s_min]#  写入文件writeTiff(MBI, geotrans, proj, MBIPath)#  原图像
filePath = r"test.jpg"
#  MBI结果
MBIPath = r"test_mbi.jpg"
#  建筑物提取结果
buildingPath = r"test_building.jpg"
#  结构元素大小最小值
s_min = 3
#  结构元素大小最大值
s_max = 20
#  测定的间隔
delta_s = 1
#  计算MBI
CalculationMBI(filePath, MBIPath, s_min, s_max, delta_s)

图片

test_mbi.jpg

MBI计算出来了以后,我们就要取阈值来提取建筑物了,阈值可以手动设置,也可以用算法自动求出阈值,这里我们采用OTSU算法[3]。

from skimage.filters import threshold_otsudef BuildingExtraction_otsu(MBIPath, buildingPath):width, height, bands, image, geotrans, proj = readTif(MBIPath)thresh = threshold_otsu(image) #返回一个阈值image[image>thresh] = 255image[image<=thresh] = 0image = image.astype(np.uint8)writeTiff(image, geotrans, proj, buildingPath)#  otsu自动计算阈值提取建筑物
BuildingExtraction_otsu(MBIPath, buildingPath)

图片

test_building.jpg

目视对照一下的话,感觉效果还不错。

参考

  1. ^Huang X and Zhang L. 2011. A multidirectional and multiscale morphological index for automatic building extraction from multispectral geoeye-1 imagery. Photogrammetric Engineering and Remote Sensing, 77(7), 721-732. [DOI: 10.14358/PERS.77.7.721]

  2. ^魏旭,高小明,岳庆兴,郭正胜.一种结合MBI和SLIC算法的遥感影像建筑物提取方法[J].测绘与空间地理信息,2019,42(10):100-103.

  3. ^otsu(大津算法)-百度百科 https://baike.baidu.com/item/otsu/16252828?fr=aladdin

来源:应用推广部

供稿:技术研发部

编辑:方梅

这篇关于python实现形态学建筑物指数MBI提取建筑物及数据获取的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/487928

相关文章

大模型研发全揭秘:客服工单数据标注的完整攻略

在人工智能(AI)领域,数据标注是模型训练过程中至关重要的一步。无论你是新手还是有经验的从业者,掌握数据标注的技术细节和常见问题的解决方案都能为你的AI项目增添不少价值。在电信运营商的客服系统中,工单数据是客户问题和解决方案的重要记录。通过对这些工单数据进行有效标注,不仅能够帮助提升客服自动化系统的智能化水平,还能优化客户服务流程,提高客户满意度。本文将详细介绍如何在电信运营商客服工单的背景下进行

基于MySQL Binlog的Elasticsearch数据同步实践

一、为什么要做 随着马蜂窝的逐渐发展,我们的业务数据越来越多,单纯使用 MySQL 已经不能满足我们的数据查询需求,例如对于商品、订单等数据的多维度检索。 使用 Elasticsearch 存储业务数据可以很好的解决我们业务中的搜索需求。而数据进行异构存储后,随之而来的就是数据同步的问题。 二、现有方法及问题 对于数据同步,我们目前的解决方案是建立数据中间表。把需要检索的业务数据,统一放到一张M

关于数据埋点,你需要了解这些基本知识

产品汪每天都在和数据打交道,你知道数据来自哪里吗? 移动app端内的用户行为数据大多来自埋点,了解一些埋点知识,能和数据分析师、技术侃大山,参与到前期的数据采集,更重要是让最终的埋点数据能为我所用,否则可怜巴巴等上几个月是常有的事。   埋点类型 根据埋点方式,可以区分为: 手动埋点半自动埋点全自动埋点 秉承“任何事物都有两面性”的道理:自动程度高的,能解决通用统计,便于统一化管理,但个性化定

python: 多模块(.py)中全局变量的导入

文章目录 global关键字可变类型和不可变类型数据的内存地址单模块(单个py文件)的全局变量示例总结 多模块(多个py文件)的全局变量from x import x导入全局变量示例 import x导入全局变量示例 总结 global关键字 global 的作用范围是模块(.py)级别: 当你在一个模块(文件)中使用 global 声明变量时,这个变量只在该模块的全局命名空

使用SecondaryNameNode恢复NameNode的数据

1)需求: NameNode进程挂了并且存储的数据也丢失了,如何恢复NameNode 此种方式恢复的数据可能存在小部分数据的丢失。 2)故障模拟 (1)kill -9 NameNode进程 [lytfly@hadoop102 current]$ kill -9 19886 (2)删除NameNode存储的数据(/opt/module/hadoop-3.1.4/data/tmp/dfs/na

异构存储(冷热数据分离)

异构存储主要解决不同的数据,存储在不同类型的硬盘中,达到最佳性能的问题。 异构存储Shell操作 (1)查看当前有哪些存储策略可以用 [lytfly@hadoop102 hadoop-3.1.4]$ hdfs storagepolicies -listPolicies (2)为指定路径(数据存储目录)设置指定的存储策略 hdfs storagepolicies -setStoragePo

Hadoop集群数据均衡之磁盘间数据均衡

生产环境,由于硬盘空间不足,往往需要增加一块硬盘。刚加载的硬盘没有数据时,可以执行磁盘数据均衡命令。(Hadoop3.x新特性) plan后面带的节点的名字必须是已经存在的,并且是需要均衡的节点。 如果节点不存在,会报如下错误: 如果节点只有一个硬盘的话,不会创建均衡计划: (1)生成均衡计划 hdfs diskbalancer -plan hadoop102 (2)执行均衡计划 hd

hdu1043(八数码问题,广搜 + hash(实现状态压缩) )

利用康拓展开将一个排列映射成一个自然数,然后就变成了普通的广搜题。 #include<iostream>#include<algorithm>#include<string>#include<stack>#include<queue>#include<map>#include<stdio.h>#include<stdlib.h>#include<ctype.h>#inclu

【C++】_list常用方法解析及模拟实现

相信自己的力量,只要对自己始终保持信心,尽自己最大努力去完成任何事,就算事情最终结果是失败了,努力了也不留遗憾。💓💓💓 目录   ✨说在前面 🍋知识点一:什么是list? •🌰1.list的定义 •🌰2.list的基本特性 •🌰3.常用接口介绍 🍋知识点二:list常用接口 •🌰1.默认成员函数 🔥构造函数(⭐) 🔥析构函数 •🌰2.list对象

【Prometheus】PromQL向量匹配实现不同标签的向量数据进行运算

✨✨ 欢迎大家来到景天科技苑✨✨ 🎈🎈 养成好习惯,先赞后看哦~🎈🎈 🏆 作者简介:景天科技苑 🏆《头衔》:大厂架构师,华为云开发者社区专家博主,阿里云开发者社区专家博主,CSDN全栈领域优质创作者,掘金优秀博主,51CTO博客专家等。 🏆《博客》:Python全栈,前后端开发,小程序开发,人工智能,js逆向,App逆向,网络系统安全,数据分析,Django,fastapi