1

I want to batch process a large number of raster files, calculate some values and properties and export these to a .csv or .xls. I managed to get a good working model that does this perfectly.

However, the next step is to export this to Python and with this I encounter some problems. I am not an expert in Python code. I also want to try this to get a little bit acquainted with it. The most important aim of exporting it is that the model should be compatible to use in different working environments.

Exporting the code is difficult because I used some tools that are not exportable to Python: Iterate Rasters, Collect Values and Calculate Values. Therefore I need to make a work-around.

What is the best way to accomplish that?

Here is the code posted on a GIS SE question.

I hope this is clear enough. Otherwise just ask for clarification.

Community
  • 1
Gerrit
  • 21
  • 3

1 Answers1

1

I managed to create a script that does the things I wanted to do. You can find the code I used below. Maybe this can be of some help to people in similar situations.

# Import necessary modules
import csv, arcpy, os, shutil
from arcpy.sa import *

# (1) -----------------------------------------------------------------------

# Set location of DEM raster files (adapt file path)
arcpy.env.workspace = "C:/blahblah/"

# Set workspace for intermediate files (adapt file path) - [CHOOSE EMPTY OR NON-EXISTENT FOLDER!]
scratchws = "C:/blahblah/"

# Set location for .csv file (adapt file path)
csvloc = "C:/blahblah/"

# (2) -----------------------------------------------------------------------

# Create scratch workspace if not present
if not os.path.exists(scratchws): os.makedirs(scratchws)
# Cleanup scratch workspace
for root, dirs, files in os.walk(scratchws):
        for a in files:
                os.unlink(os.path.join(root, a))
        for b in dirs:
                shutil.rmtree(os.path.join(root, b))

# Create .csv file location folder if not present
if not os.path.exists(csvloc): os.makedirs(csvloc)

# Enable overwrite intermediate files
arcpy.env.overwriteOutput = True

# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")
arcpy.CheckOutExtension("3D")

# Define csv file name and extension
print 'CREATE CSV FILE...'
csvfile = csvloc+'raster_properties.csv'
# Define header column names
header = ['DEM_name', 'Elevation_mean (m)', 'Elevation_min (m)', 'Elevation_max (m)',
          'Elevation_std (m)', 'Surface_area (m2)', 'Surface_area_above_mean_elevation (m2)',
          'Slope_mean (degrees)', 'Slope_min (degrees)', 'Slope_max (degrees)',
          'Slope_std (degrees)', 'Curvature_mean', 'Curvature_min', 'Curvature_max',
          'Curvature_std', 'Curvature_prof_mean', 'Curvature_prof_min', 'Curvature_prof_max',
          'Curvature_prof_std', 'Curvature_plan_mean', 'Curvature_plan_min', 'Curvature_plan_max',
          'Curvature_plan_std']
# Open csv file and write header row
print 'WRITE HEADERS...'
with open(csvfile, 'wb') as f:
        w = csv.writer(f)
        w.writerow(header)
print 'DONE'

# Create a list of GRIDs from the workspace
print 'START PROCESSING RASTERS'
rasters = arcpy.ListRasters("*", "GRID")
print 'Raster list:', rasters

# Perform all calculations for each raster (r) in the defined workspace
for r in rasters:
        # Intermediate cleanup of scratch workspace
        for root, dirs, files in os.walk(scratchws):
                for a in files:
                        os.unlink(os.path.join(root, a))
                for b in dirs:
                        shutil.rmtree(os.path.join(root, b))

        demname = str(r) # Assign raster name in string format to variable (Assign always in string format to variable from now on)
        print 'Currently processing raster:', demname

        # Calculate elevation stats and assign to variables
        arcpy.CalculateStatistics_management(r, "", "", "", "SKIP_EXISTING", "")
        elevationmean = str(arcpy.GetRasterProperties_management(r, "MEAN"))
        elevationmin = str(arcpy.GetRasterProperties_management(r, "MINIMUM"))
        elevationmax = str(arcpy.GetRasterProperties_management(r, "MAXIMUM"))
        elevationstd = str(arcpy.GetRasterProperties_management(r, "STD"))
        print 'Elevation stats (m) (mean, min, max, std):', elevationmean, elevationmin, elevationmax, elevationstd

        # Execute IsNull
        arcpy.gp.IsNull(r, scratchws+"isnull")
        # Execute SetNull
        arcpy.gp.SetNull(scratchws+"isnull", scratchws+"isnull", scratchws+"setnull", "VALUE=1")
        # Execute ZonalGeometry - calculate surface area raster
        arcpy.gp.ZonalGeometry(scratchws+"setnull", "Value", scratchws+"zonalgeomet", "AREA", r)
        arcpy.CalculateStatistics_management(scratchws+"zonalgeomet", "", "", "", "SKIP_EXISTING", "")
        # Assign surface area value to variable
        surfacearea = str(arcpy.GetRasterProperties_management(scratchws+"zonalgeomet", "MEAN"))
        print 'Surface area (m2):', surfacearea

        # Execute Zonal Statistics - calculate mean elevation raster
        arcpy.gp.ZonalStatistics(scratchws+"setnull", "Value", r, scratchws+"patch_mean", "MEAN", "DATA")
        # Execute Raster Calculator - calculate amount of cells above mean elevation
        outCon = Con((Raster(r) > Raster(scratchws+"patch_mean")), 1, 0)
        outCon.save(scratchws+"patch_m_sep")
        # Execute Set Null
        arcpy.gp.SetNull(scratchws+"patch_m_sep", scratchws+"patch_m_sep", scratchws+"above_mean", "VALUE=0")
        # Execute Zonal Geometry - calculate surface area above mean elevation
        arcpy.gp.ZonalGeometry(scratchws+"above_mean", "Value", scratchws+"area_ab_mean", "AREA", r)
        # Execute Calculate Statistics
        arcpy.CalculateStatistics_management(scratchws+"area_ab_mean", "", "", "", "SKIP_EXISTING", "")
        # Execute Get Raster Properties and assign value to variable
        surfaceareaabovemeanelevation = str(arcpy.GetRasterProperties_management(scratchws+"area_ab_mean", "MEAN"))
        # Delete patch_m_sep raster for next loop
        arcpy.Delete_management(scratchws+"patch_m_sep")
        print 'Surface area above mean elevation (m2):', surfaceareaabovemeanelevation

        # Execute Slope
        outSlope = Slope(r, "DEGREE", 1)
        # Save the output 
        outSlope.save(scratchws+"slope")
        # Calculate slope stats and assign to variables
        arcpy.CalculateStatistics_management(scratchws+"slope", "", "", "", "SKIP_EXISTING", "")
        slopemean = str(arcpy.GetRasterProperties_management(scratchws+"slope", "MEAN"))
        slopemin = str(arcpy.GetRasterProperties_management(scratchws+"slope", "MINIMUM"))
        slopemax = str(arcpy.GetRasterProperties_management(scratchws+"slope", "MAXIMUM"))
        slopestd = str(arcpy.GetRasterProperties_management(scratchws+"slope", "STD"))
        arcpy.Delete_management(scratchws+"slope")
        print 'Slope stats (degrees) (mean, min, max, std):', slopemean, slopemin, slopemax, slopestd

        # Execute Curvature (all types)
        outCurv = Curvature(r, 1, scratchws+"curv_prof", scratchws+"curv_plan")
        # Save the output 
        outCurv.save(scratchws+"curv")
        # Calculate curvature stats and assign to variables
        arcpy.CalculateStatistics_management(scratchws+"curv", "", "", "", "SKIP_EXISTING", "")
        arcpy.CalculateStatistics_management(scratchws+"curv_prof", "", "", "", "SKIP_EXISTING", "")
        arcpy.CalculateStatistics_management(scratchws+"curv_plan", "", "", "", "SKIP_EXISTING", "")
        curvmean = str(arcpy.GetRasterProperties_management(scratchws+"curv", "MEAN"))
        curvmin = str(arcpy.GetRasterProperties_management(scratchws+"curv", "MINIMUM"))
        curvmax = str(arcpy.GetRasterProperties_management(scratchws+"curv", "MAXIMUM"))
        curvstd = str(arcpy.GetRasterProperties_management(scratchws+"curv", "STD"))
        curvprofmean = str(arcpy.GetRasterProperties_management(scratchws+"curv_prof", "MEAN"))
        curvprofmin = str(arcpy.GetRasterProperties_management(scratchws+"curv_prof", "MINIMUM"))
        curvprofmax = str(arcpy.GetRasterProperties_management(scratchws+"curv_prof", "MAXIMUM"))
        curvprofstd = str(arcpy.GetRasterProperties_management(scratchws+"curv_prof", "STD"))
        curvplanmean = str(arcpy.GetRasterProperties_management(scratchws+"curv_plan", "MEAN"))
        curvplanmin = str(arcpy.GetRasterProperties_management(scratchws+"curv_plan", "MINIMUM"))
        curvplanmax = str(arcpy.GetRasterProperties_management(scratchws+"curv_plan", "MAXIMUM"))
        curvplanstd = str(arcpy.GetRasterProperties_management(scratchws+"curv_plan", "STD"))
        arcpy.Delete_management(scratchws+"curv")
        arcpy.Delete_management(scratchws+"curv_prof")
        arcpy.Delete_management(scratchws+"curv_plan")
        print 'Standard curvature stats (mean, min, max, std):', curvmean, curvmin, curvmax, curvstd
        print 'Profile curvature stats (mean, min, max, std):', curvprofmean, curvprofmin, curvprofmax, curvprofstd
        print 'Plan curvature stats (mean, min, max, std):', curvplanmean, curvplanmin, curvplanmax, curvplanstd

        # Write raster poperties to csv file
        print 'WRITING TO CSV FILE...'
        with open(csvfile, 'ab') as f:
                w = csv.writer(f)
                w.writerow([demname, elevationmean, elevationmin, elevationmax, elevationstd,
                            surfacearea, surfaceareaabovemeanelevation, slopemean, slopemin,
                            slopemax, slopestd, curvmean, curvmin, curvmax, curvstd, curvprofmean,
                            curvprofmin, curvprofmax, curvprofstd, curvplanmean, curvplanmin,
                            curvplanmax, curvplanstd])
        print 'DONE'
        print 'NEXT RASTER'
print 'FINISHED'
Gerrit
  • 21
  • 3