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

摘录于 https://blog.csdn.net/weixin_43123242/article/details/93525175

1.下载第三方包

在网址 https://www.lfd.uci.edu/~gohlke/pythonlibs/#lxml下载对应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.原图像文件&计算结果

原始图像

 

 提取边界的结果图（4领域）

 

 图3 提取边界的结果图（8领域）