How to get arcpy to access attributes of a geometry?

Question

As a novice python programmer I am struggling mightily with my code. What it is supposed to do is

1. Buffer a bunch of points with four distances (100, 200, 300, 400meters)
2. 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()

Answer 1

Create Mosaic dataset,NOT MOSAIC TO NEW RASTER. Use it as input, it will work as single raster