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
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 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) 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+".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()); showPyMessage() message = "Python Error Info: " + str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()