As a novice python programmer I am struggling mightily with my code. What it is supposed to do is
- Buffer a bunch of points with four distances (100, 200, 300, 400meters)
- Calculate zonal sum for the four distances off a raster
The problem is that the raster I need to calculate is 1m cell size for the state of California so it is actually 47 tifs that are tiled together.
I have a feature class (called CA_TileIndex) which has the tile name stored in the attribute TileName.
It works fine when the buffered zones are within the one tile. I did hack together something that used the PairWiseIntersect tool to intersect CA_TileIndex with the points but it only uses the points (not the buffered distance) and I cannot seem to figure out how the proper capture/trap when the buffered zone is in more than one tile.
buf is a geometry that I want to intersect with CA_TileIndex but I cannot use the intersect method because CA_TileIndex is a feature class, not a geometry. If I convert CA_TileIndex to geometry, I could loop through them (there are 47 tiles) but I am not clear how to get the attribute TileName from the geometry.
I am kinda stuck and am wondering if I am going about this correctly.
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
UrbanTileFC = os.path.join(WSpace,"CA_TileIndex") # This is the tile index
TSpace = r'L:\AirQuality\TreeTiles\Tile_2012' # Where the FRAP tiles reside
rasterTile = os.path.join(TSpace,"R03_C03_TreeCanopy_CA_2012.tif")
#create Tile Geometry
tileGeoms = arcpy.CopyFeatures_management(UrbanTileFC, arcpy.Geometry())
fieldsPairWise = ('TileName')
if arcpy.Exists(outFC):
arcpy.AddMessage(" Deleting " + outFC + " ...")
arcpy.Delete_management(outFC)
#pairWiseIntersect(inFC, UrbanTreeFC, outFC, fieldsPairWise)
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", "TILENAME"]
pFields=['Shape@',"D_100","D_200","D_300","D_400"]
distances=[100,200,300,400] #distances are in meters
victRaster=r'in_memory\rtem'
#victRaster = os.path.join(WSpace,"rastTemp")
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)
for g in tileGeoms:
# first check for disjoint which should make it run faster
if g.disjoint(buf):
arcpy.AddMessage("Disjoint")
else:
arcpy.AddMessage("The tile should be: ?")
#arcpy.AddMessage(" Getting tiles")
#fieldsTile= 'TileName'
#myList = sorted(unique_values(outFC,fieldsTile))
#if len(myList) == 1:
# rasterTile=os.path.join(TSpace,myList[0]+".tif")
# arcpy.AddMessage(" raster tile = " + str(rasterTile))
arcpy.Clip_management(rasterTile,envelope,victRaster,buf,"-3.402823e+038","ClippingGeometry")
extract = arcpy.RasterToNumPyArray(victRaster,"","","",-9999)
#arcpy.AddMessage(str(extract))
#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()
