Mountain vertex extraction (ArcPy implementation)

Keywords: Python arcgis Arcpy

1. Background

Mountain vertices refer to areas within a specific neighborhood analysis range where the point is higher than the surrounding points. Mountain vertices are important features of the terrain, and their distribution and density reflect the characteristics of the development of the landform, but also restrict the development of the landform. Therefore, how to extract mountain vertices correctly and effectively based on DEM data is of great importance in digital terrain analysis.

2. Purpose

By extracting and configuring contours, mountain vertices and depressions, guide readers to master the functions of contour extraction, neighborhood analysis and window calculation in spatial analysis of ArcGIS raster data, complete surface analysis of raster data.

3. Data:

1:10000 DEM data (Chp8\Ex5\) in Loess Hilly areas.

4. Requirements

(1) Using the contour extraction function in the raster data spatial analysis module, contour maps with contour distances of 15 meters and 75 meters are extracted, and contour lines are drawn according to the standard topographic map drawing method as the topographic background for mountain vertex extraction;

(2) The mountain vertices are extracted by neighborhood analysis and raster calculator.

5. Process

Load DEM data - > Load Spatial_Analyst module, Neighborhood analysis - > Map Algebra - > Reclassification - > Raster to Point

Using DEM to generate contours with intervals of 15m and 75m, a hillshade result map is generated, which forms a topographic halo map to assist in determining the location of mountain vertices. Focal statistical analysis is performed on DEM data, which is processed in a window of 21*21. After reclassifying the results with the DEM data, a raster-like Hill vertex data can be obtained. After converting the raster data into vectors, the topographic halo map is combined.Delete the unreasonable mountain vertices to get their distribution.

The flowchart is as follows:

6. Model Builder

7. ArcPy Implementation

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# 8-5 Mountain Vertex Extraction.py
# Created on: 2021-10-10 22:29:58.00000
#   (generated by ArcGIS/ModelBuilder)
# Description: 
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy
from arcpy import Raster, ListRasters, Delete_management
from arcpy.sa import Test
import os
import shutil

path = raw_input("Enter the absolute path to the folder where the data resides:").decode("utf-8")
paths = path + "\\result"
if not os.path.exists(paths):
    os.mkdir(paths)
else:
    shutil.rmtree(paths)
    os.mkdir(paths)

# Local variables:
dem = path + "\\dem"
Contour_dem15 = "Contour_dem15"
Contour_dem75 = "Contour_dem75"
HillSha_dem = "hillSha_dem"
FocalSt_dem = "FocalSt_dem"
F_dem = "F_dem"
raster_peak = "raster_peak"
vector_peak = "peak"

# Set Geoprocessing environments
print "Set Geoprocessing environments"
arcpy.env.scratchWorkspace = paths  # Set up a temporary workspace
arcpy.env.workspace = paths  # Set up workspace
arcpy.env.extent = dem  # Set Processing Range
arcpy.env.cellSize = dem  # Set Cell Size
arcpy.env.mask = dem  # Set Mask

# Process:Contour
print "Process: Contour 15 m"
arcpy.gp.Contour_sa(dem, Contour_dem15, "15", "0", "1", "CONTOUR", "")

# Process: Contour (2)
print "Process: Contour 75 m"
arcpy.gp.Contour_sa(dem, Contour_dem75, "75", "0", "1", "CONTOUR", "")

# Process:Hillshade
print "Process: Hillshade"
arcpy.gp.HillShade_sa(dem, HillSha_dem, "315", "45", "NO_SHADOWS", "1")

# Process: Focus Statistics
print "Process: Focal Statistics"
arcpy.gp.FocalStatistics_sa(dem, FocalSt_dem, "Rectangle 21 21 CELL", "MAXIMUM", "DATA", "90")

# Process:Raster Calculator
# arcpy.gp.RasterCalculator_sa("\"%FocalSt_dem%\" - \"%dem%\" == 0", F_dem)

# Process: FocalSt_dem - dem == 0
print "Process: FocalSt_dem - dem == 0"
Test(Raster(FocalSt_dem) - Raster(dem), "VALUE=0").save(F_dem)

# Process:Reclassification
print "Process: Reclassification"
arcpy.gp.Reclassify_sa(F_dem, "Value", "0 NODATA;1 1", raster_peak, "DATA")

# Process:Raster to Point
print "Process: Raster to Point"
tempEnvironment0 = arcpy.env.outputZFlag
arcpy.env.outputZFlag = "Disabled"
tempEnvironment1 = arcpy.env.outputMFlag
arcpy.env.outputMFlag = "Disabled"
arcpy.RasterToPoint_conversion(raster_peak, vector_peak, "VALUE")
arcpy.env.outputZFlag = tempEnvironment0
arcpy.env.outputMFlag = tempEnvironment1

save = ["contour_dem15", "contour_dem75", "hillsha_dem", u"peak"]
rasters = ListRasters()
for raster in rasters:
    if raster.lower() not in save:
        print u"Deleting{}layer".format(raster)
        Delete_management(raster)
print "Finished running~~~"

8. Results



Note: The number of extracted mountain vertices is limited by the size of the analysis window in block statistics. The larger the window, the fewer points will be extracted. However, if the window is too large, some important mountain vertices will be missed. For the extracted results, it can be determined manually to delete local points.

End of experiment byebye~~

Posted by dragin33 on Sun, 17 Oct 2021 09:23:15 -0700