Run Raster Geoprocessing using ArcPy#

1. Set “mesc” for ArcPy processing#

import arcpy

gdb_worksp = r"..\data\class_data.gdb"
arcpy.env.workspace = gdb_worksp

1.1 env object#

The env object exposes environment settings for geoprocessing functions in ArcGIS.

These properties can be used to retrieve (read) the current values or to set (write) them.

Geoprocessing environment settings can be thought of as additional parameters that affect a tool’s results.

1.2 Set mask#

We can set the “mask” of a geoprocessing tool using a feature class.

arcpy.env.mask = "county_boundary"

It can also be set by a raster dataset.

arcpy.env.mask = "dem"

1.3 Set extent and Read extent#

Similarly, we can set the processing “extent” using a feature class.

arcpy.env.extent = "county_boundary"

Or, a raster dataset.

arcpy.env.extent = "dem"

The extent is essentially defined by 4 numbers: top, left, bottom, right. For example, we can read the left, right, top, bottom coordinates of the extent setting.

extent_info = "Left:\t{}\nRight:\t{}\nTop:\t{}\nBottom:\t{}".format(
    arcpy.env.extent.XMin, 
    arcpy.env.extent.XMax, 
    arcpy.env.extent.YMax, 
    arcpy.env.extent.YMin
)

print(extent_info)

Left:	529453.5700000036
Right:	588493.5700000036
Top:	660647.1999999955
Bottom:	602507.1999999955

1.4 Set snap raster#

Since “snap raster” is a raster itself, so it can only be set by a raster dataset.

arcpy.env.snapRaster = "dem"

1.5 Set cell size#

We can set the output raster’s cell size by another raster.

arcpy.env.cellSize = "dem"

Or a number. Note that this number is in the default linear unit defined by the raster’s spatial reference.

arcpy.env.cellSize = 90

2. Run Raster Functions using ArcPy#

# setup "mesc"

cntbnd_fc = "county_boundary"
dem_rast = "dem"

arcpy.env.mask = cntbnd_fc
arcpy.env.extent = cntbnd_fc
arcpy.env.snapRaster = None
arcpy.env.cellSize = 30

2.1 Setup output geodatabase#

Recall the ways we can define an output feature class here.

import os

output_gdb = r"..\data\output_data.gdb"
os.path.join(output_gdb, 'out_rast')

'..\\data\\output_data.gdb\\out_rast'

output_gdb + "\\" + 'out_rast' # with string concatnation

'..\\data\\output_data.gdb\\out_rast'

"{}\{}".format(output_gdb, 'out_rast') # with .format function

'..\\data\\output_data.gdb\\out_rast'

output_gdb = os.path.abspath(output_gdb)

2.2 Euclidean Distance#

?arcpy.sa.EucDistance

school_fc = "schools"

school_dist_rast = arcpy.sa.EucDistance(school_fc)   # without setting cell_size
school_dist_rast.save(os.path.join(output_gdb, 'school_dist_rast'))

2.3 Slice#

?arcpy.sa.Slice

arcpy.env.snapRaster = school_dist_rast

school_dist_1_9 = arcpy.sa.Slice(
    school_dist_rast,
    number_zones=9, 
    slice_type="NATURAL_BREAKS",
    base_output_zone=1
)

school_dist_1_9.save(output_gdb + "\\" + "school_dist_1_9")

2.4 Reclassify#

  • RemapValue object takes a list of lists, with the inner lists being composed of 2 parts

    • oldValue

    • newValue

  • RemapRange object takes a list of lists, with the inner lists being composed of 3 parts

    • startValue

    • endValue

    • newvalue

?arcpy.sa.Reclassify

dist_remap = arcpy.sa.RemapValue(
    [[1, 9], [2, 8], [3, 7], [4, 6],
     [6, 4], [7, 3], [8, 2], [9, 1]]
)

school_dist_9_1 = arcpy.sa.Reclassify(school_dist_1_9, "Value", dist_remap)
school_dist_9_1.save("{}\{}".format(output_gdb, "school_dist_9_1"))

2.5 Map Algebra directly operated on rasters#

school_dist_wgted = school_dist_9_1*0.3 + school_dist_1_9*0.7
school_dist_wgted.save(os.path.join(output_gdb, "school_dist_wgted"))

?arcpy.sa.Int

# without defining a variable to store the output raster

arcpy.sa.Int(school_dist_wgted).save(
    os.path.join(output_gdb, "school_dist_wgted_int")
)