1

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[0]
            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[0] for r in arcpy.da.SearchCursor(UrbanAreaFC, ['SHAPE@'])][0]
                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()[2])[0]; showPyMessage()
    message = "Python Error Info: " +  str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()

closed as off-topic by PolyGeo Dec 17 '15 at 2:29

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "Questions seeking help to debug/write/improve code must include the desired behavior, a specific problem or error and the shortest code necessary to reproduce it in the question itself. Providing a clear problem statement and evidence of a code attempt will help others to help you. See: How to create a Minimal, Complete, and Verifiable example." – PolyGeo
If this question can be reworded to fit the rules in the help center, please edit the question.

  • 2
    The entire script would need to be re-written. This is too broad of a question. – Richard Law Dec 17 '15 at 2:29
  • As a first step to focussing your question I recommend that you decide whether you really want to attempt a re-write in pyQGIS, or whether fixing it in ArcPy may be a better way forward. – PolyGeo Dec 17 '15 at 2:31
  • Hmm. Well fixing it in ArcPy is not really an option. As I stated I am not at their location so I cannot really diagnose why the script is failing. – user918967 Dec 17 '15 at 5:18