June2025/KD_REM at main · JVanDerHout/June2025 · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114

import arcpy
from arcpy.sa import *
import os


# Check out the Spatial Analyst extension
arcpy.CheckOutExtension("Spatial")

# Set environment settings
arcpy.env.overwriteOutput = True

# ---------------- User-Defined Parameters ---------------- #

# Path to the DEM raster (LiDAR DEM)
dem_path = r"\\neo\Terra\ClientFiles\Q-T\TuusiWanaDesign_220220\GIS\Elevation_Data\ProposedSurface_240923_9cfsEG_burn.tif"

# Path to the channel line feature class (CL within the geodatabase)
channel_line = r"\\neo\Terra\ClientFiles\Q-T\TuusiWanaDesign_220220\GIS\PRO\TuusiWana_CottonWoodZones\TuusiWana_CottonWoodZones.gdb\CL"

# Output workspace geodatabase
output_workspace = r"\\neo\Terra\ClientFiles\Q-T\TuusiWanaDesign_220220\GIS\PRO\TuusiWana_CottonWoodZones\TuusiWana_CottonWoodZones.gdb"

# Point spacing distance along the channel line (e.g., approximate channel width)
point_spacing = 30  # in map units (e.g., meters)

# Search radius for the Kernel Density tool (should cover the floodplain area of interest)
search_radius = 330  # in map units (e.g., meters)

# --------------------------------------------------------- #

# Set the workspace environment to the output geodatabase
arcpy.env.workspace = output_workspace

# Ensure the output workspace exists
if not arcpy.Exists(output_workspace):
    arcpy.CreateFileGDB_management(os.path.dirname(output_workspace), os.path.basename(output_workspace))

# Set the cell size to match the DEM cell size
dem_raster = arcpy.Raster(dem_path)
cell_size = dem_raster.meanCellWidth  # Assuming square cells

# Set environment cell size to match DEM
arcpy.env.cellSize = cell_size

# ------------------- Processing Steps -------------------- #

# Step KD-Step2b: Generate points along the channel line
print("Generating points along the channel line...")
channel_points = os.path.join(output_workspace, "channel_points")

arcpy.management.GeneratePointsAlongLines(
    Input_Features=channel_line,
    Output_Feature_Class=channel_points,
    Point_Placement="DISTANCE",
    Distance=f"{point_spacing} Feet",
    Include_End_Points=True
)

# Step KD-Step2c: Extract elevation values to channel points
print("Extracting elevation values to channel points...")
arcpy.sa.ExtractMultiValuesToPoints(
    in_point_features=channel_points,
    in_rasters=[[dem_path, "ELEV"]],
    bilinear_interpolate_values="NONE"
)

# Step KD-Step3: Create the point density raster using Kernel Density
print("Creating the point density raster...")
kd_point_density = KernelDensity(
    in_features=channel_points,
    population_field="NONE",
    cell_size=cell_size,
    search_radius=search_radius,
    area_unit_scale_factor="SQUARE_MAP_UNITS",
    out_cell_values="DENSITIES",
    method="PLANAR"
)

point_density_raster = os.path.join(output_workspace, "point_density")
kd_point_density.save(point_density_raster)

# Step KD-Step4: Create the stream elevation density raster
print("Creating the stream elevation density raster...")
kd_stream_elev_density = KernelDensity(
    in_features=channel_points,
    population_field="ELEV",
    cell_size=cell_size,
    search_radius=search_radius,
    area_unit_scale_factor="SQUARE_MAP_UNITS",
    out_cell_values="DENSITIES",
    method="PLANAR"
)

stream_elev_density_raster = os.path.join(output_workspace, "stream_elev_density")
kd_stream_elev_density.save(stream_elev_density_raster)

# Step KD-Step5: Create the detrended DEM by dividing the two rasters
print("Creating the detrended DEM...")
detrended_dem = arcpy.sa.Divide(kd_stream_elev_density, kd_point_density)
detrended_dem_raster = os.path.join(output_workspace, "detrended_dem")
detrended_dem.save(detrended_dem_raster)

# Step KD-Step6: Generate the Relative Elevation Model (REM)
print("Generating the Relative Elevation Model (REM)...")
raw_dem = arcpy.Raster(dem_path)
rem = arcpy.sa.Minus(raw_dem, detrended_dem)
rem_raster = os.path.join(output_workspace, "kd_rem_PRLowFlow_330r")
rem.save(rem_raster)

# Check in the Spatial Analyst extension
arcpy.CheckInExtension("Spatial")

print("Processing complete. The REM has been created at:", rem_raster)