基于Python和GDAL实现栅格数据中相邻地物边界提取

欢迎客官移步到简研–简简单单科研中学习和交流

对输入的影像计算指定像元值四邻域内指定像元值的像元数并输出成新文件

此代码块是打包成exe的源文件

1、导入相关python包

在网址 https://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal 下载对应python版本的whl文件。如，GDAL‑3.0.0‑cp38‑cp38m‑win32.whl

pip install numpy
pip install GDAL‑3.0.0‑cp38‑cp38m‑win32.whl

2、代码

from osgeo import gdal, gdalconst
from osgeo.gdalconst import *
import numpy as np
import sys
np.set_printoptions(threshold = 1e6) 
#设置输入和输出参数
#参数1：输入待计算边界的原始TIFF影像
#参数2：计算的结果影像文件（TIFF文件）
#参数3：参数1中的待计算的焦点像元值
#参数4：用于计算边界像元的邻域算子窗口大小
#参数5：参数1中的待计算的邻域像元值
Input = sys.argv[1].replace('\\','/')
FocusPoint = int(sys.argv[2])
ZoneNear = int(sys.argv[3])
NeiOpe = int(sys.argv[4])
Output = sys.argv[5].replace('\\','/')
#读取栅格数据
ds = gdal.Open(Input,GA_ReadOnly)
if ds is None:
    print(Input)
cols = ds.RasterXSize
rows = ds.RasterYSize
geotransform = ds.GetGeoTransform()
geoProjection = ds.GetProjection()
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]
band = ds.GetRasterBand(1)
data = band.ReadAsArray(0, 0, cols, rows)
data = data.astype(np.int)
Ordata = np.array(data,dtype = int)
#基于原始数据构造二元值
UniqueValue = np.unique(Ordata)#计算唯一像元值
OnlyFocusPoint = np.where(Ordata == FocusPoint, 0, -1)
OnlyZoneNear = np.where(Ordata == ZoneNear, 2, 0)
FZ = OnlyFocusPoint + OnlyZoneNear
ReData = np.where(FZ == -1, 0, FZ)
#拼接数据
row0 = np.zeros([1,cols], dtype = int)
col0 = np.zeros([rows+2,1], dtype = int)
rowPinRow = np.r_[row0,ReData,row0]
rowPinCol = np.c_[col0,rowPinRow,col0]
DataPin = rowPinCol
rowsPin = np.shape(DataPin)[0]
colsPin = np.shape(DataPin)[1]
outData = np.zeros([rowsPin,colsPin],dtype = np.int)
#构造切片
if NeiOpe == 8: 
    #8邻域,不包括中心像元
    outData[1:rowsPin-1,1:colsPin-1] = (DataPin[0:rowsPin-2,0:colsPin-2] + DataPin[0:rowsPin-2,1:colsPin-1] + DataPin[0:rowsPin-2,2:colsPin] + DataPin[1:rowsPin-1,0:colsPin-2] + DataPin[1:rowsPin-1,2:colsPin] + DataPin[2:rowsPin,0:colsPin-2] + DataPin[2:rowsPin,1:colsPin-1] + DataPin[2:rowsPin,2:colsPin])
elif NeiOpe == 4:
    #4邻域，不包括中心像元
    outData[1:rowsPin-1,1:colsPin-1] = (DataPin[0:rowsPin-2,1:colsPin-1] + DataPin[1:rowsPin-1, 0:colsPin-2] + DataPin[1:rowsPin-1,2:colsPin] + DataPin[2:rowsPin,1:colsPin-1])
else:
    print('Only 4 or 8')
ResultData = outData[1:rowsPin-1,1:colsPin-1]
#构造淹没
Mask = np.where(Ordata == FocusPoint, 0, np.nan)
EdgeData = np.array(Mask + ResultData)
#新建栅格用于存放EdgeData
driver = gdal.GetDriverByName("GTiff")
outDataset = driver.Create(Output, cols, rows, 1, gdal.GDT_Int16)
outDataset.SetGeoTransform(geotransform) 
outDataset.SetProjection(geoProjection)
outBand = outDataset.GetRasterBand(1)
outBand.WriteArray(EdgeData)
outBand.SetNoDataValue(0)
outDataset.FlushCache()
#至此计算指定像元值的焦点像元邻域中特地像元值的像元个数计算完成
#若计算出具体的边界长度，可用pixelWidth或pixelHeight乘以EdgeData计算即可
print('Done')

3、打包成exe可执行文件

使用PyInstaller可以将Python程序打包成Windows可执行程序，此处不写详细步骤，参考：https://www.cnblogs.com/yu2000/p/3797654.html

4、原图像文件&计算结果

图1 提取边界的原图

图2 提取边界的结果图（4领域）
注释：像元值为1，表示图1中的建设用地4领域内共有1个耕地像元。

图3 提取边界的结果图（8领域）
注释：像元值为1，表示图1中的建设用地8领域内共有1个耕地像元。

这里分界线

安利个看美剧背单词，练习听说的工具|美剧词典

美剧词典是一款短情景形式的英语学习工具。首先呢，颜值是十分高的，界面简洁，配色清新；其次，有才华啊，目前共收录了近10万条字幕数据，10万条短情景和80万条词汇。短情景均为几秒的无字幕原声经典美剧视频片段，不占用多少时间就能让你get到原汁原味的单词含义和用处啦。美剧词典现包括**“查询”、“单词本”、“音频"和"我的”**四项功能。