Find out the growth area of a precious medicinal material (ArcPy Implementation)

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)