Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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.

share|improve this question

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 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

# 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',
# Open csv file and write header row
print 'WRITE HEADERS...'
with open(csvfile, 'wb') as f:
        w = csv.writer(f)
print 'DONE'

# Create a list of GRIDs from the workspace
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, scratchws+"isnull")
        # Execute SetNull"isnull", scratchws+"isnull", scratchws+"setnull", "VALUE=1")
        # Execute ZonalGeometry - calculate surface area raster"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"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)"patch_m_sep")
        # Execute Set Null"patch_m_sep", scratchws+"patch_m_sep", scratchws+"above_mean", "VALUE=0")
        # Execute Zonal Geometry - calculate surface area above mean elevation"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
        print 'Surface area above mean elevation (m2):', surfaceareaabovemeanelevation

        # Execute Slope
        outSlope = Slope(r, "DEGREE", 1)
        # Save the output"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"))
        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"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"))
        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'
share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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