0

I'm using ArcGIS Pro and I need to calculate the proportion of pixels within 500m for a specified value. I have a raster with coded values of either 0 or 1. For each pixel, I want to calculate the proportion of 1s within a 500m radius. The focal statistics function seems to be the closest tool to accomplish this, but I'm not getting the results I expect.

# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *

# Set local variables
inRaster = "Hazard_Rating"
neighborhood = NbrCircle(500, "MAP")  # Define a circular neighborhood with 500-meter radius

# Define the constant circle neighborhood with a radius of 500 meters (7854 pixels)
circle_radius = 500
circle_area = 7854
circle_neighborhood = NbrCircle(circle_radius, "CELL")


# Execute FocalStatistics to Calculate the Sum of ones
outFocalStatistics = FocalStatistics(inRaster, circle_neighborhood, "SUM")

# Calculate the proportion of 1s within the fixed circle neighborhood
proportionRaster = outFocalStatistics / circle_area

This seems like it should work but I'm getting values above 1 in proportionRaster, and when I check outFocalStatistics the raster attribute table is empty.

4
  • a) use mean as suggested by @user2856 or recalculate circle area
    – FelixIP
    Commented Aug 12, 2023 at 4:42
  • Can you explain this a little more? why would calculating the mean be the same as the percentage? how will recalculating the circle area help? thank you!
    – seak23
    Commented Aug 12, 2023 at 8:22
  • Mean IS what you call percentage. If you insist on using total, recalculate area to number of cells when using map units = 715398/cellSize/cellSize
    – FelixIP
    Commented Aug 12, 2023 at 21:44
  • @seak23 did this answer your question, if so gis.stackexchange.com/help/someone-answers. If not, how could the answer be improved?
    – user2856
    Commented Mar 28 at 2:30

1 Answer 1

1

Your current calculation will only be correct if your cell size is 1x1m.

Either:

  • Use FocalStatistics(inRaster, circle_neighborhood, "MEAN")

    Because you have a grid of 0 and 1 values, mean = sum/count = proportion of 1's in neighbourhood.

    By using mean, you also let ArcGIS handle the edge cases where the focal neighbourhood is not a complete circle. E.g the output at any of the corners will be the proportion of 1s in a quarter circle.

  • If you want to treat edge/corner pixels the same (i.e assume the circle can extend outside the bounds of the raster), calculate the circle area in grid cells/pixels so you are calculating the proportion using the same units

    from arcpy.sa import (FocalStatistics, Raster)
    from math import pi
    
    inRaster = Raster("Hazard_Rating")
    
    circle_radius = 500  # m
    circle_area = pi * circle_radius**2
    cell_area = inRaster.meanCellWidth * inRaster.meanCellHeight
    # Estimate the number of grid cells in the circle. This is just approximate
    # See https://pro.arcgis.com/en/pro-app/latest/arcpy/spatial-analyst/nbrcircle-class.htm
    circle_count = circle_area // cell_area  # // is integer division and truncates the result
    
    # Define a circular neighborhood with 500m radius because we specfied "MAP" not "CELL"
    circle_neighborhood = NbrCircle(circle_radius, "MAP")  
    
    # Execute FocalStatistics to Calculate the Sum of ones
    outFocalStatistics = FocalStatistics(inRaster, circle_neighborhood, "SUM")
    
    # Calculate the proportion of 1s within the fixed circle neighborhood
    proportionRaster = outFocalStatistics / circle_count
    
    

Note that the above estimate circle_count won't give exactly the number of pixels in circle_neighborhood. There are ways to manually calculate this, however the simplest way to do this is just to expand the analysis extent by +- radius, fill the outside NoData areas with zeros and crop back to the original extent after calculating the focal mean, e.g.

  from arcpy.sa import (FocalStatistics, Raster)
  from math import pi
  
  inRaster = Raster("raster_one.tif")
  
  radius = 500  # m
  
  # Define a circular neighborhood with 500m radius because we specfied "MAP" not "CELL"
  circle_neighborhood = NbrCircle(radius, "MAP")  
  
  extent = inRaster.extent
  
  xmin, ymin, xmax, ymax = extent.XMin - radius, extent.YMin - radius, extent.XMax + radius, extent.YMax + radius
  
  with arcpy.EnvManager(extent=arcpy.Extent(xmin, ymin, xmax, ymax)):
  
      # Execute FocalStatistics to Calculate the Sum of ones
      outFocalStatistics = FocalStatistics(Con(IsNull(inRaster), 0, inRaster), circle_neighborhood, "MEAN")
  
  with arcpy.EnvManager(extent=extent):
      outRaster = ApplyEnvironment(outFocalStatistics)
 

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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