利用代码实现山脊线、山谷线的提取（arcpy版）

山脊线和山谷线的实现代码其实大部分都是一样的，下面我会在适当的地方指出来,具体效果根据需要再改改吧

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)