I have a land use layer with a map of the United States overlaid on top. I'd like to determine the percentage of different land uses in each state. Would the best way be to extract by mask by using each state and then calculating within each state?
I am very much a novice GIS user so any assistance/guidance would be incredibly useful! I am using ArcGIS 10.
If you have a Java background then python should be easy.. here's something that I put together:
import sys, os, arcpy
InRaster = sys.argv[1]
InShape = sys.argv[2]
SplitField = sys.argv[3]
OutFolder = sys.argv[4]
# simple switches based on the output, change as needed
# note that python is case sensitive (true != True) so in
# the case of booleans False is a constant as is True
# but true, TRUE or any other combination is not
do_zoning = False # do Zonal Statistics as Table
do_extract = True # do Extract by Mask to create extracted rasters
# check out spatial analyst if you're running from a DOS prompt
if arcpy.CheckExtension("Spatial") == "Available":
arcpy.AddMessage( "Checking out Spatial")
arcpy.CheckOutExtension("Spatial")
else:
arcpy.AddMessage( "Unable to get spatial analyst extension")
sys.exit(0)
F = arcpy.ListFields(InShape,SplitField)
if not F:
arcpy.AddError("Field not found " + SplitField)
sys.exit(-1)
if not os.path.exists(OutFolder):
arcpy.AddError("Output path does not exist")
sys.exit(-2)
Values = list() # new empty list
# step through the features adding each unique value
# to the list such that we can then step through later
arcpy.AddMessage("Obtaining unique values")
with arcpy.da.SearchCursor(InShape,[ SplitField ]) as SCur:
for SRow in SCur:
if not SRow[0] in Values:
Values.append(SRow[0])
for ThisValue in Values:
# make a feature layer based on just this value
defQ = "%s = %s" % (SplitField,ThisValue)
arcpy.MakeFeatureLayer_management(InShape,"Layer",defQ)
# set output(s)
outDBF = "%s\\ZonalStat_%s.dbf" % (OutFolder,ThisValue)
outRaster = "%s\\Extract_%s.tif" % (OutFolder,ThisValue)
if do_zoning:
arcpy.AddMessage("Zoning " + ThisValue)
arcpy.sa.ZonalStatisticsAsTable("Layer",SplitField,InRaster,outDBF)
if do_extract:
# extract the raster using the layer to an arcpy Raster object
py_raster = arcpy.sa.ExtractByMask(InRaster,"Layer")
py_raster.save(outRaster) # save the raster object to a raster
# clean up
arcpy.Delete_management("Layer") # just the layer, not the data, so this name can be reused
arcpy.CheckInExtension("Spatial") # release spatial analyst extension (if running from DOS)
Have a read of Adding a Script Tool if you want to run it from catalog, otherwise you can run from DOS (cmd) window.
As it is written the code will do an extraction but not a zonal statistics, although the code is there I've turned it off with a boolean variable (do_zoning).
