導航

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

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')