As a novice python programmer I am struggling mightily with my code. What it is supposed to do is
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 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 polygonsUrbanTileFC = os.path.join(WSpace,"CA_TileIndex") # This is the tile indexTSpace = r'L:\AirQuality\TreeTiles\Tile_2012' # Where the FRAP tiles residerasterTile = os.path.join(TSpace,"R03_C03_TreeCanopy_CA_2012.tif")#create Tile GeometrytileGeoms = 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 metersvictRaster=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) # 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()
أكثر...
- Buffer a bunch of points with four distances (100, 200, 300, 400meters)
- Calculate zonal sum for the four distances off a raster
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 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 polygonsUrbanTileFC = os.path.join(WSpace,"CA_TileIndex") # This is the tile indexTSpace = r'L:\AirQuality\TreeTiles\Tile_2012' # Where the FRAP tiles residerasterTile = os.path.join(TSpace,"R03_C03_TreeCanopy_CA_2012.tif")#create Tile GeometrytileGeoms = 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 metersvictRaster=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) # 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()
أكثر...