找出某名珍貴藥材的生長區域(ArcPy實現)
阿新 • • 發佈:2021-10-12
某種珍貴藥材生長於山區,通過研究瞭解到這種藥材生長具有嚴格的生長條件(見下文)。為了能更好地保護該藥材的生長環境,現在需要使用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 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
Collection:因為這些東西是非常簡單的。不要抱怨自己學不會,那是因為你沒有足夠用心。