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, timeimport itertools as itfrom arcpy import envenv.overwriteoutput=Trueenv.workspace=r'in_memory'WSpace = r'C:\Users\Don\Documents\AirQuality\AirBufferAnalysis.gdb' # Where the file geodatabase is with your points in itinFC = os.path.join(WSpace, "TBTestPoint_2") # Point file with the geocoded addressesoutFC = os.path.join(WSpace,"TB_Processed") # Final Output NameUrbanAreaFC = os.path.join(WSpace,"CA_UrbanAreas") # This is the urban census polygonsrasterTile = 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 metersvictRaster=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) # 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()
أكثر...
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, timeimport itertools as itfrom arcpy import envenv.overwriteoutput=Trueenv.workspace=r'in_memory'WSpace = r'C:\Users\Don\Documents\AirQuality\AirBufferAnalysis.gdb' # Where the file geodatabase is with your points in itinFC = os.path.join(WSpace, "TBTestPoint_2") # Point file with the geocoded addressesoutFC = os.path.join(WSpace,"TB_Processed") # Final Output NameUrbanAreaFC = os.path.join(WSpace,"CA_UrbanAreas") # This is the urban census polygonsrasterTile = 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 metersvictRaster=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) # 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()
أكثر...