0

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()
1

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

  • I tried that already. The problem was that the mosaic tries to make a 3-band image and the line 'reduced=filter(lambda x: x not in [-9999],chain)' bombs – user918967 May 13 '15 at 19:59
  • Yes, I have a MOSAIC DATASET. Running my code gives me a <type 'exceptions.ValueError'> The truth value of an array with more than one element is ambiguous. Use a.any() or a.all(). I think it is using the single raster, but I think with three bands the RasterToNumPyArray returns a larger array? – user918967 May 13 '15 at 20:39
  • script designed for 1 band raster.period – FelixIP May 13 '15 at 20:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.