Sometime ago, I asked a question (Arcpy to append results of pivot table) with the result being python 2.7 script that works in ArcGIS 10.3. It is fairly simple - it iterates over all the points in a feature class and then for 100, 200, 300, 400 meters around the point performs a focal sum on a raster mosaic dataset (that is only 1s and 0s). I took the answer and modified it only slightly (to perform focal sum rather than focal mean).
Unfortunately, I have a colleague that is using ArcGIS 10.1 and this script does not work on his system. It could be the backwards compatibility issue (10.1 vs 10.3) and it also could be the raster data mosaic (which is terribly finicky to set up).
Either way, since I am not onsite in their office I figured that rather than pound my head on a wall since his shop uses QGIS (which I have never used), the python code needs to be tweaked and updated to methods QGIS uses.
Having not used QGIS in the past, I thought I could turn to the group for assistance. What changes to the existing code are necessary?
import arcpy, traceback, os, sys, time import itertools as it from arcpy import env env.overwriteoutput=True env.workspace=r'in_memory' WSpace = r'C:\Users\Don\Documents\AirQuality\AirBufferAnalysis.gdb' # Where the file geodatabase is with your points in it inFC = os.path.join(WSpace, "TBTestPoint_2") # Point file with the geocoded addresses outFC = os.path.join(WSpace,"TB_Processed") # Final Output Name UrbanAreaFC = os.path.join(WSpace,"CA_UrbanAreas") # This is the urban census polygons rasterTile = os.path.join(WSpace,"TreeData") if arcpy.Exists(outFC): arcpy.AddMessage(" Deleting " + outFC + " ...") arcpy.Delete_management(outFC) arcpy.CopyFeatures_management(inFC, outFC) arcpy.AddField_management(outFC, "D_100", "DOUBLE") arcpy.AddField_management(outFC, "D_200", "DOUBLE") arcpy.AddField_management(outFC, "D_300", "DOUBLE") arcpy.AddField_management(outFC, "D_400", "DOUBLE") pFields=['Shape@',"D_100","D_200","D_300","D_400"] distances=[100,200,300,400] #distances are in meters victRaster=r'in_memory\rtem' try: def showPyMessage(): arcpy.AddMessage(str(time.ctime()) + " - " + message) pgon=arcpy.Geometry() with arcpy.da.UpdateCursor(outFC,pFields) as cursor: m=0 for row in cursor: point=row arcpy.AddMessage("Starting Point Number: " + str(m)) for i in range(4): buf=point.buffer(distances[i]) # make sure the buffer is totally within Urban Areas urban_area_geom = [r for r in arcpy.da.SearchCursor(UrbanAreaFC, ['SHAPE@'])] if buf.within(urban_area_geom): anExtent=buf.extent envelope='%f %f %f %f' %(anExtent.XMin, anExtent.YMin, anExtent.XMax, anExtent.YMax,) if arcpy.Exists(victRaster): arcpy.Delete_management(victRaster) arcpy.Clip_management(rasterTile,envelope,victRaster,buf,"-3.402823e+038","ClippingGeometry") extract = arcpy.RasterToNumPyArray(victRaster,"","","",-9999) extract = arcpy.RasterToNumPyArray(victRaster,"","","",0) chain = it.chain.from_iterable(extract) reduced=filter(lambda x: x not in [-9999],chain) mySum=sum(reduced) #The output is the square meters of trees in the buffered zone row[1+i]=mySum else: row[1+i]=-9999 m+=1 cursor.updateRow(row) #arcpy.AddMessage("Finished Point Number: " + str(m)) except: message = "\n*** PYTHON ERRORS *** "; showPyMessage() message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()); showPyMessage() message = "Python Error Info: " + str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()