<commit_before>### Calculate the treeline from a DEM and a forest raster
### Output is a raster with points along the treeline
### calculates for each pixel if it a) the it is forest b) the highest forest pixel in the vicinity
### and c) if there exists non forest land at higher altitude in the vicinity


import os, sys, arcpy

from arcpy import env
from arcpy.sa import *


##################################


# Set Environments
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")
arcpy.env.workspace = 'D:/Test/Michael/Wald_Jan/Scratch/'

# Inputs / Outputs
inDGM = 'D:/Test/Michael/Wald_Jan/DGM_Mongolei_SRTM_90.tif'
inForest = 'D:/Test/Michael/Wald_Jan/Wald_Raster_90m.tif'
outFile = 'D:/Test/Michael/Wald_Jan/Waldgrenzen/Waldgrenze.tif'

# Vicinity to check
forestRadius = 5 # is there forest at higher altitude in this radius
nonForest = 5    # is there a higher elevation at all in this radius

# Check vicinity for altitude differences, if diff is below threshold, terrain is assumed flat and point excluded
flatRadius = 5  #cells
flatDiff = 100    #meters

# Altitude difference in meter, mini
alt_diff = 50


##################################

dgm = Raster(inDGM)
forest = Raster(inForest)
forest1 = forest+1 # change forest raster values from 0 to 1

# Get altitude values only where forest exisits
alt_forst = SetNull(forest, dgm, "VALUE <> 0")

# Create Raster Layers from input data
arcpy.MakeRasterLayer_management(inDGM, "dgm_layer")
arcpy.MakeRasterLayer_management(alt_forst, "alt_forest_layer")
arcpy.MakeRasterLayer_management(forest1, "forest1_layer")

# create Focal Statistics, each pixel in the output raster gets the maximum value of the cell in radius
alt_forstFoc = FocalStatistics("alt_forest_layer", NbrAnnulus(0, forestRadius, "CELL"), 
                               "MAXIMUM", "DATA")
dgmFoc = FocalStatistics("dgm_layer", NbrAnnulus(0, nonForest, "CELL"), 
                               "MAXIMUM", "DATA")

# Focal Stats to decrease points
altDiffFoc = FocalStatistics("dgm_layer", NbrAnnulus(0, nonForest, "CELL"), 
                               "RANGE", "DATA")
minForestFoc = FocalStatistics("forest1_layer", NbrAnnulus(0, 10, "CELL"), 
                               "SUM", "DATA")

# create new raster with 1 where highest altitude in forest in vicinity and where higher elev. exisits
outCon = Con( (alt_forst == alt_forstFoc) &
              (dgmFoc-alt_diff > dgm) &
              (altDiffFoc > 50) &
              (minForestFoc > 10), 1)

outCon.save(outFile[:-4] + "_" + str(forestRadius) + "_" +
            str(nonForest) + "_" + str(alt_diff) + "_" +
            str(flatRadius) + "_" + str(flatDiff) + outFile[-4:])


# Delete files created on HDD
arcpy.Delete_management(alt_forstFoc)
arcpy.Delete_management(dgmFoc)
arcpy.Delete_management(altDiffFoc)


arcpy.CheckInExtension("Spatial")<commit_msg>lower tree line added, slope exclusion added, cleaned<commit_after>### Calculate the treeline from a DEM and a forest raster
### Output is a raster with points along the treeline
### calculates for each pixel if it a) the it is forest b) the highest forest pixel in the vicinity
### and c) if there exists non forest land at higher altitude in the vicinity


import os, sys, arcpy

from arcpy import env
from arcpy.sa import *


##################################
 # ENVIRONMENTS AND PARAMETERS #
##################################


# Set Environments
arcpy.env.overwriteOutput = True
#arcpy.env.cellSize = 'MINOF'
#arcpy.env.snapRaster = 'D:/Test/Michael/Wald_Jan/Wald_Raster_90m.tif'
arcpy.env.workspace = 'D:/Test/Michael/Wald_Jan/Scratch/'

arcpy.CheckOutExtension("Spatial")



# Inputs / Outputs
inDGM = 'D:/Test/Michael/Wald_Jan/DGM_Mongolei_SRTM_90.tif' # DEM file in meters altitude
inForest = 'D:/Test/Michael/Wald_Jan/Wald_Raster_90m.tif'   # Forest file, 0 = Forest, NoData = no forest
inSlope = 'D:/Test/Michael/Wald_Jan/Slope_WGS_degree.tif'       # Slope in Degree per Pixel

outFile = 'D:/Test/Michael/Wald_Jan/Waldgrenzen/Waldgrenze.tif' # Basic output name, will be extended by parameters


# Vicinity to check
forestRadius = 5 # is there forest at higher altitude in this radius
nonForest = 5    # is there a higher elevation at all in this radius

# Altitude difference in meter, minimum slope
alt_diff = 50
minSlope = 2


##################################
 # MAIN CALCULATION PART #
##################################


dgm = Raster(inDGM)
slope = Raster(inSlope)
forest = Raster(inForest)
forest1 = forest+1 # change forest raster values from 0 to 1


# Get altitude values only where forest exisits
alt_forst = SetNull(forest, dgm, "VALUE <> 0")


# Create Raster Layers from input data
arcpy.MakeRasterLayer_management(inDGM, "dgm_layer")
arcpy.MakeRasterLayer_management(alt_forst, "alt_forest_layer")


# create Focal Statistics, each pixel in the output raster gets the max/min value of the cell in radius
alt_max_forstFoc = FocalStatistics("alt_forest_layer", NbrAnnulus(0, forestRadius, "CELL"), 
                               "MAXIMUM", "DATA")
alt_min_forstFoc = FocalStatistics("alt_forest_layer", NbrAnnulus(0, forestRadius, "CELL"), 
                               "MINIMUM", "DATA")
dgmFoc_max = FocalStatistics("dgm_layer", NbrAnnulus(0, nonForest, "CELL"), 
                               "MAXIMUM", "DATA")
dgmFoc_min = FocalStatistics("dgm_layer", NbrAnnulus(0, nonForest, "CELL"), 
                               "MINIMUM", "DATA")


# create new raster with 1 where highest altitude in forest in vicinity and where higher elev. exisits
outCon = Con( (alt_forst == alt_max_forstFoc) &
              (dgmFoc_max - alt_diff > dgm) &
              (slope > minSlope), 1, # upper tree boundary 

              Con( (alt_forst == alt_max_forstFoc) &
                   (dgmFoc_min+alt_diff > dgm) &
                   (slope > minSlope), 2) # lower tree boundary
              
              )


##################################
 # OUTPUT CREATION AND DELETION #
##################################

outCon.save(outFile[:-4] + "_" + str(forestRadius) + "_" +
            str(nonForest) + "_" +
            str(minSlope) + outFile[-4:])


# Delete files created on HDD
arcpy.Delete_management(alt_max_forstFoc)
arcpy.Delete_management(alt_min_forstFoc)
arcpy.Delete_management(dgmFoc_max)
arcpy.Delete_management(dgmFoc_min)


arcpy.CheckInExtension("Spatial")
{'fraction_non_alphanumeric': 0.13711324014449505, 'fraction_numerical': 0.009423590387937804, 'mean_word_length': 3.7275426874536004, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3650745', 'n_tokens_mistral': 2200, 'n_tokens_neox': 2033, 'n_words': 570}