1

I have a forested area, that suffered different disturbances over multiple years. Each year, disturbance data was collected on different forest extend and some disturbances can overlap (beetle infestation or clear-cut).

My goal is to clean and summarize the dataset, i.e. unify the extent and allow the forest area to be disturbed only once over my time period. Simply, I wish to clip the disturbance data to the available forest area in current year, and update the extend of forest cover after every disturbance.

I am trying to apply arcpy to automate the process.

Pseudocode:

  • read forest area & list yearly disturbances
  • create for loop to loop through list of disturbances
    • clip disturbance by available forest
    • update available forest by erasing the disturbance from current forest cover
    • do the same for next disturbance

I thought this could be fairy simple to do. However, my script either return empty dataset, erase the initial forest or says "that output and input data equals".

My desired output is the set of cleaned (not overlapping disturbances limited on forest) disturbances. I am trying to automate the process, as I have 3 different disturbances types over multiple years, together 18 disturbances featureclasses.

How to update the forest polygon and spatially restrict every disturbance?

import arcpy
import os

# Allow files to overwrite
arcpy.env.overwriteOutput = True

# Define directories
path = "U:\\project\\analyzed"
arcpy.env.workspace = os.path.join(path, "analyzed.gdb")
outWD = os.path.join(path, "cleaned.gdb")
aoi = "aoi_fin_03"

# Read all disturbance data
disturbLs = arcpy.ListFeatureClasses()

# Intersect available forest with disturbances in each year,
# update forest cover 
for dist in disturbLs:
    # Clip the disturbance to forest extend
    outClip = os.path.join(outWD, dist)
    arcpy.Clip_analysis(dist, aoi, outClip)

    # Update available forest cover by erasing current disturbance
    updAoi = "in_memory\layer" + aoi
    arcpy.CopyFeatures_management(aoi, updAoi)
    arcpy.Erase_analysis(updAoi, outClip, aoi)

#Delete in memory files to avoid the lock
arcpy.Delete_management("in_memory")
  • How many different disturbance datasets do you have? – BERA May 10 at 12:22
  • I have beetle infestation, clear cut and windthorow, varying in years. In total 18 feature classes, which I order using this script: disturbLs.sort(key=lambda x:[int(x.split('_')[1]),x.split('_')[0]]) (not shown in example to keep it simple) – maycca May 10 at 12:26
1

This might be oversimplifying the problem but a single Union and Select could work:

  1. Union stands and all disturbances
  2. Select where FID_stands>-1 and all disturbances FID<0
  3. Export selection

enter image description here

  • thank you for your answer! However, my forest does not consist of stands. It is derived from the raster forest cover data and as such I have only 1 forest polygon. That's why I need to do clips or intersections, and erasing. Do you think that the logic of my pseudocode is correct? – maycca May 10 at 13:14
  • It should work with only one forest polygon, does not matter – BERA May 10 at 13:30
  • Hi BERA, I have tried your approach. It seems interesting, but pretty complicated: I have 18 disturbances in total, which makes to manually select 18 fields in a query. Also, as an optput I expect to get unique clean disturbances, again 18 in total. Those I would like to later group by year, and calculate sum of yearly non overlapping disturbances. – maycca May 10 at 14:14
  • You can construct the Query using arcpy and listfields , no need to do it manually – BERA May 13 at 7:08
0

Actually, my problem was that for some features clip tool generated empty output. This lead to to further generate empty forest class etc.

So here, I check if the count of rows in a feature class and than run my loop only on those features:

if arcpy.management.GetCount(dist)[0] > "0":
    print dist

The whole code:

import arcpy
import os

# Allow files to overwrite
arcpy.env.overwriteOutput = True

# Define directories
path = "U:\\projects\\analyzed"
arcpy.env.workspace = os.path.join(path, "trySubset.gdb")
outWD = os.path.join(path, "cleanedSub.gdb")
aoi = "aoi_fin_03"


# Read all disturbance data
disturbLs = arcpy.ListFeatureClasses()

# For loop
for dist in disturbLs:

    # Check if the feature class is not empty
    if arcpy.management.GetCount(dist)[0] > "0":

        # clip disturbance to forest area, 
        # save the output
        outClipDist = os.path.join(outWD, dist)
        arcpy.Clip_analysis(dist, aoi, outClipDist)

        #  save the updated forest to verify it
        updAOI = os.path.join(outWD, "aoi_" + dist)

        # Erase the current disturbance from the available forest
        arcpy.Erase_analysis(aoi, outClipDist, updAOI)

        # update the forest cover to use at the beginning of for loop
        aoi = updAOI  

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.