1

I'm trying to write a code in Python for several rasters, however I'm wondering if there is a faster way to do this.

I'm generating several rasters from a conditional statement for each value in the raster file. Then from these rasters I'm doing focal statistics. I'm saving all the interim rasters from the conditional statement and using them in focal statistics. Is there a way to avoid that? Can I make the code below more automatic?

# Import arcpy module
import arcpy
from arcpy.sa import *
#
arcpy.CheckOutExtension("spatial")
arcpy.env.overwriteOutput = True
arcpy.env.workspace = "D:\\_Work\\Test-Model.gdb"
#
CropMap = "RasterFile"
#
#
# Step 1
# Select crops
select_crop = "select_crop_136"
arcpy.gp.Con_sa(CropMap, "1", select_crop, "0", "Value = 136")
#
select_crop = "select_crop_146"
arcpy.gp.Con_sa(CropMap, "1", select_crop, "0", "Value = 146")
#
select_crop = "select_crop_147"
arcpy.gp.Con_sa(CropMap, "1", select_crop, "0", "Value = 147")
#
#
# Total Crop
tot_crop = "tot_crop_136"
arcpy.gp.Con_sa(CropMap, "1", tot_crop, "1", "Value = 136")
#
tot_crop = "tot_crop_146"
arcpy.gp.Con_sa(CropMap, "1", tot_crop, "1", "Value = 146")
#
tot_crop = "tot_crop_147"
arcpy.gp.Con_sa(CropMap, "1", tot_crop, "1", "Value = 147")
#
#
# Step 2
# Focal Statistics
#
neighborhood = NbrRectangle(1000, 1000, "MAP")
#
inRaster_sel_136 = "select_crop_136"
foc_select_136 = FocalStatistics(inRaster_sel_136, neighborhood, "SUM","")
#
inRaster_tot_136 = "tot_crop_136"
foc_tot_136 = FocalStatistics(inRaster_tot_136, neighborhood, "SUM", "")
#
# and so forth for all the rasters in Step 1
3
  • 1
    Why aren't you using the new (at V10) Raster objects and arcpy.sa.Con, in step 1 like in Step 2? It's not immediately clear what you're trying to achieve in your result (outFocalStatistics) but either way you're not saving it, just overwriting the existing raster in the step 2 process so there's no real output generated that I can see. Dec 20 '17 at 1:23
  • In addition to @MichaelStimson comment, your tot_crop rasters are all the same as long as true == false_true, which is a constant raster for your raster extent.
    – fatih_dur
    Dec 20 '17 at 1:25
  • Thank you Michael and faith_dur for your answers. I edited the code according to your comments: (i) replaced the outFocalStatistics with a unique name and (ii) replaced "1" and "0" in the code
    – egi007
    Jan 4 '18 at 20:05
0

This is untested, but around about how I would tackle this. I am not considering your logic as mentioned by others, just setting it up for you. Also, make it common practice to never have special characters in your file names that aren't underscores (e.g. '-' in your workspace).

There is room for improvement on this, but this should work for you.

import arcpy
from arcpy.sa import *

arcpy.CheckOutExtension("spatial")
arcpy.env.overwriteOutput = True
arcpy.env.workspace = "D:\\_Work\\Test-Model.gdb"

CropMap = "RasterFile"

crop = [136, 146, 147]
neighborhood = NbrRectangle(1000, 1000, "MAP")

for e in crop:
    select_crop = Con(Raster(crop_map) == e, 1, 0)
    total_crop = Con(Raster(crop_map) == e, 1, 0)
    arcpy.gp.FocalStatistics_sa(select_crop, "select_{}".format(str(e)), neighborhood , "SUM", "DATA")
    arcpy.gp.FocalStatistics_sa(total_crop, "total_{}".format(str(e)), neighborhood , "SUM", "DATA")
3
  • Thanks so much @David Campbell! The Con part runs fine. I get an error at the FocalStatistics part. Maybe there is a problem with the output raster name: "select_{}".format(str(e))?
    – egi007
    Jan 4 '18 at 19:10
  • what does the error say exactly?
    – D_C
    Jan 4 '18 at 22:13
  • 1
    Now it works with FocalStatistics(select_crop, neighborhood, "SUM", "") Thanks
    – egi007
    Jan 10 '18 at 2:00

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.