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

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

import os
import gdal
import arcpy
from arcpy import env
import numpy as npdef read_img(filename):dataset=gdal.Open(filename)im_width = dataset.RasterXSizeim_height = dataset.RasterYSizeim_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_datadef write_img(filename, im_proj, im_geotrans, im_data):if 'int8' in im_data.dtype.name:datatype = gdal.GDT_Byteelif 'int16' in im_data.dtype.name:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float32if len(im_data.shape) == 3:im_bands, im_height, im_width = im_data.shapeelse: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 datasetdef 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_datadem_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.0correct_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] = 0correct_data[correct_data > 60.0] = 1temp1 = os.path.join(temp_path, 'temp1.tif')write_img(temp1, im_proj, im_geotrans, correct_data)dem_mean_c[dem_mean_c >= 0] = 0dem_mean_c[dem_mean_c < 0] = 1temp2 = os.path.join(temp_path, 'temp2.tif')write_img(temp2, im_proj, im_geotrans, correct_data)arr = correct_data + dem_mean_carr[arr < 2] = 0arr[arr >= 2] = 1write_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 = Trueenv.workspace = temp_pathsg_extract(dem_path, valley_out, temp_path)

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

利用水文分析提取山脊线山谷线
1 流程图利用水文分析提取山脊线及山谷线,山脊线相当于分水线,山谷线相当于山谷线.分水线是水流的起源点,这些栅格的水流方向只存在流出方向而不存在流入方向,所以汇流累积量为零.通过对零值的提取就可以得 ...
基于matlab山脊线,教你如何利用水文，分析提取山脊线山谷线
原标题:教你如何利用水文,分析提取山脊线山谷线利用水文分析提取山脊线及山谷线,山脊线相当于分水线,山谷线相当于山谷线.分水线是水流的起源点,这些栅格的水流方向只存在流出方向而不存在流入方向,所以汇流 ...
matlab 山脊提取,ArcGIS中利用水文分析提取山脊线山谷线
1 流程图利用水文分析提取山脊线及山谷线,山脊线相当于分水线,山谷线相当于山谷线.分水线是水流的起源点,这些栅格的水流方向只存在流出方向而不存在流入方向,所以汇流累积量为零.通过对零值的提取就可以得 ...
基于matlab山脊线,山脊线山谷线提取实验报告.doc
山脊线山谷线提取实验报告实验内容描述: 山脊线和山谷线构成了地形起伏变化的分界线(骨架线),因此它对于地形地貌研究具有重要意义:另一方面,对于水文物理过程研究而言,由于山脊.山谷分别代表示分水性与汇 ...
Arcgis中山脊线，山谷线的提取，以及流域的分割
目录实验数据下载:https://pan.baidu.com/s/1jfTv5LXBDaabguJsPTFJEw 提取码:4g1b Task 1:利用Ex1中的数据,练习提取不同位置的地形剖面线: ...
ArcGIS山脊线、山谷线和山顶点的提取(附练习数据下载)
特征地形要素,主要是指对地形在地表的空间分布特征具有控制作用的点.线或面状要素.特征地形要素构成地表地形与起伏变化的基本框架.与地形指标的提取主要采用小范围的邻域分析不同的是,特征地形要素的提取更多地 ...
利用水文分析方法提取山脊、山谷线
1.背景作为地形特征线的山音线.山谷线对地形.地貌具有一定的控制作用.它们与山顶点.谷底点以及鞍部点等一起构成了地形起伏变化的骨架结构.同时由于山登线具有分水性,山谷线具有合水性特征,使得它们在地形 ...
【ArcGIS风暴】水文分析模块实验：山脊线和山谷线提取
实验平台:ArcGIS 9.3 实验目的:学习和掌握山脊线和山谷线提取的原理及方法实验要求:利用ArcGIS水文分析模块提取样区的山脊线和山谷线实验数据:Ex1 实验步骤: 1.正负地形的提取 ( ...
[空间分析] DEM提取山脊线与山谷线
文章目录基于图像处理技术的原理基于地形表面几何形态分析原理的算法基于地形表面流水物理模拟分析原理的算法基于地形表面几何形态分析和流水物理模拟分析相结合平面曲率与坡位组合法水文分析法 [山脊 ...

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

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

最新文章

热门文章