1, Background
A kind of precious medicinal material grows in mountainous areas. Through research, it is known that this kind of medicinal material has strict growth conditions. In order to better protect the growth environment of the medicinal material, it is necessary to use GIS spatial analysis method to find out the suitable growth area of the medicinal material, so as to provide basis for the protection of the species.
2, Data
(1) Mountain contour data contour. shp;
(2) Annual average temperature and annual total precipitation data collected from observation points in mountainous area climate. txt.
3, Growth conditions of medicinal materials
Please determine the area suitable for planting this kind of medicinal material according to the following conditions, and make a thematic map to give the suitable planting area.
(1) This kind of medicinal material generally grows in the area close to both sides of the valley (generally no more than 500m);
(2) This medicine likes Yang;
(3) The growth climate environment is that the annual average temperature is 10 ° ~ 12 °;
(4) The total annual precipitation is 550 ~ 680mm.
4, Process
DEM data is generated from the mountain contour data. Hydrological analysis is carried out based on DEM to extract gully network (catchment threshold can be set according to needs); slope direction data is extracted based on DEM to reclassify Yin and Yang slopes.
The annual average temperature and annual total precipitation data collected from observation points are used for surface interpolation respectively to generate annual average temperature grid data and annual total precipitation grid data (the interpolation method is selected according to needs). The area with annual average temperature of 10 ℃ and 12 ℃ and the area with annual total precipitation of 550680mm are extracted.
The area meeting the above four conditions is comprehensively superimposed and analyzed to obtain the area suitable for the growth of the medicinal material, and the thematic map is made to calculate the area of the suitable area.
The process is shown in the figure below:
###5, Model builder
6, ArcPy implementation
Note: when creating tin and converting it into grid dem, I don't know why it always reports spatial reference parameter errors. The same parameter settings can be done in ArcMap and model builder.
At present, no solution has been found. Little partners who know how to solve it are welcome to leave a message in the comment area~
Therefore, when running the following code, you need to create tin with ArcMap or model builder, convert it into grid dem, and copy the dem to the working path.
# -*- coding: utf-8 -*- # --------------------------------------------------------------------------- # 13-4 find out the generation area of a precious medicinal material.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("Please enter the absolute path of the folder where the data is located:").decode("utf-8") # Start timing 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" # The length of the grid base name cannot exceed 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: create TIN # print "Process: create TIN" # arcpy.CreateTin_3d(tin, "Unknown", "{}.shp CONTOUR Soft_Line CONTOUR".format(contour), "DELAUNAY") # Process: TIN to grid # print "Process: TIN to grid" # arcpy.TinRaster_3d(tin, dem, "FLOAT", "LINEAR", "CELLSIZE 100", "1") arcpy.env.extent = dem arcpy.env.cellSize = dem arcpy.env.mask = dem # Process: filling depression print "Process: Filling depression" arcpy.gp.Fill_sa(dem, Fill_dem, "") # Process: flow direction print "Process: flow to" arcpy.gp.FlowDirection_sa(Fill_dem, FlowDir_Fill, "NORMAL", Descent_data, "D8") # Process: aspect print "Process: Slope direction" arcpy.gp.Aspect_sa(dem, Aspect_dem, "PLANAR", "METER") # Process: reclassification print "Process: Reclassification" arcpy.gp.Reclassify_sa(Aspect_dem, "value", "90 270 1", sunny_slope, "NODATA") # Process: creating XY event layers print "Process: establish XY Event Layer " 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: inverse distance weight method print "Process: Inverse distance weight method" arcpy.gp.Idw_sa(climate_Layer, "temperature", temperature, "100", "2", "VARIABLE 12", "") # Process: grid calculator print "Process: Grid calculator" # arcpy.gp.RasterCalculator_sa("(\"%temperature%\" >= 10) & (\"%temperature%\" <= 12)", proper_temp) ((Raster(temperature) >= 10) & (Raster(temperature) <= 12)).save(proper_temp) # Process: inverse distance weight method (2) print "Process: Inverse distance weight method (2)" arcpy.gp.Idw_sa(climate_Layer, "rainfall", precipitation, "100", "2", "VARIABLE 12", "") # Process: grid calculator (2) print "Process: Grid calculator (2)" # arcpy.gp.RasterCalculator_sa("(\"%precipitation%\" >= 550) & (\"%precipitation%\" <= 680)", proper_prec) ((Raster(precipitation) >= 550) & (Raster(precipitation) <= 680)).save(proper_prec) # Process: Traffic print "Process: flow" arcpy.gp.FlowAccumulation_sa(FlowDir_Fill, FlowAcc_Flow, "", "FLOAT", "D8") # Process: greater than or equal to print "Process: Greater than or equal to" arcpy.gp.GreaterThanEqual_sa(FlowAcc_Flow, constant1, stream) # Process: reclassification (2) print "Process: Reclassification (2)" arcpy.gp.Reclassify_sa(stream, "VALUE", "0 NODATA;1 1", Reclass_stre1, "DATA") # Process: Euclidean distance print "Process: Euclidean distance " arcpy.gp.EucDistance_sa(Reclass_stre1, distance, "", dem, Direction_data, "PLANAR", "", Inverse_data) # Process: less than or equal to print "Process: Less than or equal to" arcpy.gp.LessThanEqual_sa(distance, constant2, buffer) # Process: grid calculator (3) print "Process: Grid calculator (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"Deleting{}layer ".format(raster) arcpy.Delete_management(raster) # End timing time_end = time.time() # Calculation time time_all = time_end - time_start print time.asctime() print "Execution complete!>>><<< Total time{:.0f}branch{:.2f}second".format(time_all // 60, time_all % 60)
7, Results
This experiment, too tmex o(╥﹏╥╥) O