山脊线和山谷线的实现代码其实大部分都是一样的,下面我会在适当的地方指出来,具体效果根据需要再改改吧
import os
import gdal
import arcpy
from arcpy import env
import numpy as np
def read_img(filename):
dataset=gdal.Open(filename)
im_width = dataset.RasterXSize
im_height = dataset.RasterYSize
im_geotrans = dataset.GetGeoTransform()
im_proj = dataset.GetProjection()
im_data = dataset.ReadAsArray(0,0,im_width,im_height)
del dataset
return im_proj,im_geotrans,im_width, im_height,im_data
def write_img(filename, im_proj, im_geotrans, im_data):
if 'int8' in im_data.dtype.name:
datatype = gdal.GDT_Byte
elif 'int16' in im_data.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
if len(im_data.shape) == 3:
im_bands, im_height, im_width = im_data.shape
else:
im_bands, (im_height, im_width) = 1,im_data.shape
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)
dataset.SetGeoTransform(im_geotrans)
dataset.SetProjection(im_proj)
if im_bands == 1:
dataset.GetRasterBand(1).WriteArray(im_data)
else:
for i in range(im_bands):
dataset.GetRasterBand(i+1).WriteArray(im_data[i])
del dataset
def get_reverse(img_path,reverse_path):
im_proj,im_geotrans,im_width, im_height,im_data = read_img(img_path)
dem_reverse = im_data.max() - im_data
dem_reverse = dem_reverse.astype(np.float32)
write_img(reverse_path, im_proj, im_geotrans, dem_reverse)
def correct_slope_c(slope,slope_reverse,correct_slope_out):
im_proj,im_geotrans,im_width, im_height,slope = read_img(slope)
im_proj,im_geotrans,im_width, im_height,slope_reverse = read_img(slope_reverse)
correct_slope = ((slope + slope_reverse) - abs(slope - slope_reverse)) / 2.0
correct_slope = correct_slope.astype(np.float32)
write_img(correct_slope_out, im_proj, im_geotrans, correct_slope)
def sg_extract(dem_ras, valley_out, temp_path):
aspect_temp = os.path.join(temp_path, 'aspect.tif')
arcpy.Aspect_3d(in_raster=dem_ras, out_raster=aspect_temp)
slope_temp = os.path.join(temp_path, 'slope.tif')
arcpy.Slope_3d(in_raster=aspect_temp,out_raster=slope_temp,output_measurement="DEGREE",z_factor="1")
dem_reverse = os.path.join(temp_path, 'dem_reverse.tif')
get_reverse(dem_ras, dem_reverse)
aspect_reverse_temp = os.path.join(temp_path, 'aspect_reverse.tif')
arcpy.Aspect_3d(in_raster=dem_reverse, out_raster=aspect_reverse_temp)
slope_reverse_temp = os.path.join(temp_path, 'slope_reverse.tif')
arcpy.Slope_3d(in_raster=aspect_reverse_temp,out_raster=slope_reverse_temp,output_measurement="DEGREE",z_factor="1")
correct_slope_out = os.path.join(temp_path, 'correct_slope_temp.tif')
correct_slope_data = correct_slope_c(slope_temp,slope_reverse_temp,correct_slope_out)
dem_mean = os.path.join(temp_path, 'dem_mean.tif')
arcpy.gp.FocalStatistics_sa(dem_ras, dem_mean, "Rectangle 25 25 CELL","MEAN","DATA")
im_proj,im_geotrans,im_width, im_height,dem_data = read_img(dem_ras)
im_proj,im_geotrans,im_width, im_height,dem_mean_data = read_img(dem_mean)
dem_mean_c = dem_data - dem_mean_data
#(公式("correct_data ">70)&("dem_mean_c ">0))和山谷线(公式("correct_data ">70)&("dem_mean_c "<0)
#correct_data代表坡向变率,70是要自己设定的阈值,可以更改,dem_mean_c 是正负地形分布区域,是通过邻域分析的焦点统计得到的,具体可以自己去搜下arcgis怎么操作
im_proj,im_geotrans,im_width, im_height, correct_data = read_img(correct_slope_out)
correct_data[correct_data <= 60.0] = 0
correct_data[correct_data > 60.0] = 1
temp1 = os.path.join(temp_path, 'temp1.tif')
write_img(temp1, im_proj, im_geotrans, correct_data)
dem_mean_c[dem_mean_c >= 0] = 0
dem_mean_c[dem_mean_c < 0] = 1
temp2 = os.path.join(temp_path, 'temp2.tif')
write_img(temp2, im_proj, im_geotrans, correct_data)
arr = correct_data + dem_mean_c
arr[arr < 2] = 0
arr[arr >= 2] = 1
write_img(valley_out, im_proj, im_geotrans, arr)
if __name__ == "__main__":
dem_path = 'xxx/t_dem.tif'
valley_out = 'xxx/t_dem_re.tif'
temp_path = 'xxx/temp_dem/'
arcpy.CheckOutExtension('Spatial')
arcpy.env.overwriteOutput = True
env.workspace = temp_path
sg_extract(dem_path, valley_out, temp_path)