找出某名珍貴葯材的生長區域（ArcPy實現）

一、背景

某種珍貴葯材生長於山區，通過研究了解到這種葯材生長具有嚴格的生長條件。為了能更好地保護該葯材的生長環境，現在需要使用GIS空間分析方法，將葯材適合生長區域找出來，以便為該物種保護提供依據。

二、數據

(1) 山區等高線數據contour. shp;

(2)山區觀測點采集的年平均溫度和年總降水數據climate. txt。

三、葯材生長條件

請依據以下條件，確定此區域適合種植這種葯材的范圍，並制作專題圖，給出適宜種植的面積。

(1)這種葯材一般生長在溝谷兩側較近的區域（一般不超過500m);

(2)這種葯材喜陽;

(3)生長氣候環境為年平均溫度為10°～12°;

(4)年總降水量為550~～680mm。

四、流程

利用該山區等高線數據生成DEM數據。基於DEM進行水文分析，提取溝谷網絡(匯水閾值可根據需要自行設定);基於DEM提取坡向數據，重分類划分陰陽坡。

利用觀測點采集的年平均溫度和年總降水數據分別進行表面內插，生成年平均溫度柵格數據和年總降水柵格數據（插值方法根據需要自行選擇)。提取年平均溫10℃_{12℃的區域和年總降水為550}680mm的區域。

綜合疊加分析滿足上述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 module
import arcpy
from arcpy.sa import Raster
import os
import shutil
import time

print 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 environments
print "Set Geoprocessing environments"
arcpy.env.scratchWorkspace = paths
arcpy.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 = dem
arcpy.env.cellSize = dem
arcpy.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_start
print time.asctime()
print "執行完畢！>>><<< 共耗時{:.0f}分{:.2f}秒".format(time_all // 60, time_all % 60)

七、結果

這節實驗，太tmex o(╥﹏╥)o