找出某名珍贵药材的生长区域（ArcPy实现）

找出某名珍贵药材的生长区域（ArcPy实现）

某种珍贵药材生长于山区，通过研究了解到这种药材生长具有严格的生长条件（见下文)。为了能更好地保护该药材的生长环境，现在需要使用GIS空间分析方法，将药材适合生长区域找出来，以便为该物种保护提供依据。

一、背景

某种珍贵药材生长于山区，通过研究了解到这种药材生长具有严格的生长条件。为了能更好地保护该药材的生长环境，现在需要使用GIS空间分析方法，将药材适合生长区域找出来，以便为该物种保护提供依据。

二、数据

(1) 山区等高线数据contour. shp;

(2)山区观测点采集的年平均温度和年总降水数据climate. txt。

三、药材生长条件

请依据以下条件，确定此区域适合种植这种药材的范围，并制作专题图，给出适宜种植的面积。

(1)这种药材一般生长在沟谷两侧较近的区域（一般不超过500m);

(2)这种药材喜阳;

(3)生长气候环境为年平均温度为10°～12°;

(4)年总降水量为550~～680mm。

四、流程

利用该山区等高线数据生成DEM数据。基于DEM进行水文分析，提取沟谷网络(汇水阈值可根据需要自行设定);基于DEM提取坡向数据，重分类划分阴阳坡。

利用观测点采集的年平均温度和年总降水数据分别进行表面内插，生成年平均温度栅格数据和年总降水栅格数据（插值方法根据需要自行选择)。提取年平均温10℃12℃的区域和年总降水为550680mm的区域。

综合叠加分析满足上述4个条件的区域，得到适合该药材生长的区域，并制作专题图,计算该适合区域的面积。

流程如下图所示：

五、模型构建器

六、ArcPy实现

注：创建tin并转成栅格dem，不知道为什么老是报空间参考参数错误。一样的参数设置，在ArcMap和模型构建器又可以做得出来。 目前还未找到解决办法，知道如何解决的小伙伴，欢迎评论区留言哈~

所以，在运行下面代码时，需要先用ArcMap或模型构建器创建tin并转成栅格dem，把dem复制到工作路径。

# -*- coding: utf-8 -*-# ---------------------------------------------------------------------------# 13-4 找出某种珍贵药材的生成区域.py# Created on: 2021-10-12 20:20:20.00000# (generated by ArcGIS/ModelBuilder)# Description: # ---------------------------------------------------------------------------# Import arcpy moduleimport arcpyfrom arcpy.sa import Rasterimport osimport shutilimport timeprint time.asctime()path = raw_input("请输入数据所在文件夹的绝对路径：").decode("utf-8")# 开始计时time_start = time.time()paths = path + "\\result"if not os.path.exists(paths): os.mkdir(paths)else: shutil.rmtree(paths) os.mkdir(paths)# Local variables:# contour = path + "\\contour"climate_txt = path + "\\climate.txt"tin = "tin"dem = path + "\\dem"constant1 = "200"constant2 = "500"Fill_dem = "Fill_dem"Descent_data = "Descent_data"Aspect_dem = "Aspect_dem"sunny_slope = "sunny_slope"climate_Layer = "climate_Layer"temperature = "temperature"proper_temp = "proper_temp"precipitation = "precipitate" # 格网基本名称的长度不能超过了 13。proper_prec = "proper_prec"FlowDir_Fill = "FlowDir_Fill"FlowAcc_Flow = "FlowAcc_Flow"stream = "stream"Reclass_stre1 = "Reclass_stre1"distance = "distance"buffer = "buffer"proper_area = "proper_area"Direction_data = "Direct_data"Inverse_data = "Inverse_data"# Set Geoprocessing environmentsprint "Set Geoprocessing environments"arcpy.env.scratchWorkspace = pathsarcpy.env.workspace = paths# Process: 创建 TIN# print "Process: 创建 TIN"# arcpy.CreateTin_3d(tin, "Unknown", "{}.shp CONTOUR Soft_Line CONTOUR".format(contour), "DELAUNAY")# Process: TIN 转栅格# print "Process: TIN 转栅格"# arcpy.TinRaster_3d(tin, dem, "FLOAT", "LINEAR", "CELLSIZE 100", "1")arcpy.env.extent = demarcpy.env.cellSize = demarcpy.env.mask = dem# Process: 填洼print "Process: 填洼"arcpy.gp.Fill_sa(dem, Fill_dem, "")# Process: 流向print "Process: 流向"arcpy.gp.FlowDirection_sa(Fill_dem, FlowDir_Fill, "NORMAL", Descent_data, "D8")# Process: 坡向print "Process: 坡向"arcpy.gp.Aspect_sa(dem, Aspect_dem, "PLANAR", "METER")# Process: 重分类print "Process: 重分类"arcpy.gp.Reclassify_sa(Aspect_dem, "value", "90 270 1", sunny_slope, "NODATA")# Process: 创建 XY 事件图层print "Process: 创建 XY 事件图层"arcpy.MakeXYEventLayer_management(climate_txt, "x", "y", climate_Layer, "{B286C06B-0879-11D2-AACA-00C04FA33C20};281094.800114045 3873927.75678257 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision", "")# Process: 反距离权重法print "Process: 反距离权重法"arcpy.gp.Idw_sa(climate_Layer, "温度", temperature, "100", "2", "VARIABLE 12", "")# Process: 栅格计算器print "Process: 栅格计算器"# arcpy.gp.RasterCalculator_sa("(\"%temperature%\" >= 10) & (\"%temperature%\" <= 12)", proper_temp)((Raster(temperature) >= 10) & (Raster(temperature) <= 12)).save(proper_temp)# Process: 反距离权重法 (2)print "Process: 反距离权重法 (2)"arcpy.gp.Idw_sa(climate_Layer, "降雨", precipitation, "100", "2", "VARIABLE 12", "")# Process: 栅格计算器 (2)print "Process: 栅格计算器 (2)"# arcpy.gp.RasterCalculator_sa("(\"%precipitation%\" >= 550) & (\"%precipitation%\" <= 680)", proper_prec)((Raster(precipitation) >= 550) & (Raster(precipitation) <= 680)).save(proper_prec)# Process: 流量print "Process: 流量"arcpy.gp.FlowAccumulation_sa(FlowDir_Fill, FlowAcc_Flow, "", "FLOAT", "D8")# Process: 大于等于print "Process: 大于等于"arcpy.gp.GreaterThanEqual_sa(FlowAcc_Flow, constant1, stream)# Process: 重分类 (2)print "Process: 重分类 (2)"arcpy.gp.Reclassify_sa(stream, "VALUE", "0 NODATA;1 1", Reclass_stre1, "DATA")# Process: 欧氏距离print "Process: 欧氏距离"arcpy.gp.EucDistance_sa(Reclass_stre1, distance, "", dem, Direction_data, "PLANAR", "", Inverse_data)# Process: 小于等于print "Process: 小于等于"arcpy.gp.LessThanEqual_sa(distance, constant2, buffer)# Process: 栅格计算器 (3)print "Process: 栅格计算器 (3)"# arcpy.gp.RasterCalculator_sa("\"%sunny_slope%\" * \"%proper_temp%\" * \"%proper_prec%\"* \"%buffer%\"", proper_area)(Raster(sunny_slope) * Raster(proper_temp) * Raster(proper_prec) * Raster(buffer)).save(proper_area)save = ["tin", "dem", "proper_area"]rasters = arcpy.ListRasters()for raster in rasters: if raster.lower() not in save: print u"正在删除{}图层".format(raster) arcpy.Delete_management(raster)# 结束计时time_end = time.time()# 计算所用时间time_all = time_end - time_startprint time.asctime()print "执行完毕！>>><<< 共耗时{:.0f}分{:.2f}秒".format(time_all // 60, time_all % 60)

七、结果

这节实验，太tmex o(╥﹏╥)o

箴言：因为这些东西是非常简单的。不要抱怨自己学不会，那是因为你没有足够用心。

版权声明：本文内容由网络用户投稿，版权归原作者所有，本站不拥有其著作权，亦不承担相应法律责任。如果您发现本站中有涉嫌抄袭或描述失实的内容，请联系我们jiasou666@gmail.com 处理，核实后本网站将在24小时内删除侵权内容。