I'm writing an Arcpy program to process data in a layer of several hundred points ('ptsLyr', from file 'ptsFile'), each of which has its own ID ('RecordID') and an associated value, 'radius', that I want to use to create a distance buffer around that point as a proxy for a trade area; at the location of the point, the value of the raster will be 1, while at the edge of the radius and beyond it will be 0, with a linear decrease from center to edge. Ultimately, I want to sum all of the individual buffers to get a distance-weighted measure of how many points contribute to each cell in the raster. Because of the number of points I don't want to create hundreds of individual raster files.
As an example, consider a cell that is located between two points, A and B, such that it is 4000 m from A and 2000 m from B, which have radii of 8000 m and 10000 m, respectively. The contribution to the cell from point A will be (8000-4000)/8000 = 0.5, as the cell is 1/2 way between point A and its perimeter. The contribution from point B will be (10000-2000)/10000 = 0.8, since it is 1/5 of the way from B's center to its perimeter. If the cell is 6000 m from another point, C, that has a radius of 3000 m, C will not contribute to the calculation as any value outside 3000 m from C will be 0. At the cell, then, the 'cumulative' value of the overlapping buffers is 0.5 + 0.8 = 1.3.
I have to iterate through the points individually because the distance function EucDistance() calculates the distance from each cell to the closest point, which is not what I want. Also, I apparently have to create a new layer from each point as I process it because EucDistance() doesn't seem to respect selections of individual points made by SelectLayerByAttribute_management(). EucDistance() will set values outside the radius to 'No Data' but I can easily enough reclassify those to zero once I have the final, summed raster.
The map coordinates are projected (NAD83 UTM Zone 17N) and all values are in metres.
import os import arcpy import arcpy.mapping as amap from arcpy.sa import * mapFile = r"C:\Project Files\project.mxd" ptsFile = r"C:\Project Files\data\points.shp" sumFile = r"C:\Project Files\data\Project.gdb\sumbuffer" tmpFile = r"C:\Project Files\data\one_point.shp" mxd = amap.MapDocument(mapFile) df = amap.ListDataFrames(mxd) mlayers = amap.ListLayers(mxd) for lyr in mlayers: if lyr.name == os.path.split(ptsFile)[:-4]: ptsLyr = lyr arcpy.env.extent = arcpy.Extent(500000.0, 4705300.0, 775850.0, 4988910.0) counter = 0 # counter variable ptCursor = arcpy.SearchCursor(ptsLyr, "", "", "RecordID; radius") # get records for pt in ptCursor: # iterate through each point in ptsLyr itemID = pt.getValue("RecordID") radius = pt.getValue("radius") # define query to select each point by RecordID query = "\"RecordID\" = " + str(itemID) # create new layer with single selected point for further processing arcpy.FeatureClassToFeatureClass_conversion(ptsLyr.name, os.path.split(tmpFile), os.path.split(tmpFile)[:-4], query) # calculate distance to point out to value of radius distR = EucDistance(tmpFile, radius, 50) # cell size 50 m dist1 = radius - distR # invert values to give 0 at edge dist2 = dist1 / radius # rescale values to give 1 at center # add each raster to one that stores all previous output if counter == 0: sumRas = dist2 else: sumRas = dist2 + sumRas counter += 1 # increment counter variable sumRas.save(sumFile) # save the final summed raster in gdb
At the moment I have the basics working but am getting errors associated with either: the raster layers/raster math. I suspect the problems have something to do with my environment settings, raster objects vs layers, files vs geodatabases, and/or the (im)mutability of rasters.
If I do the raster math in the format 'dist1 = radius - distR', as above, I get the error
RuntimeError: ERROR 010240: Could not save raster dataset to C:\Users\JohnR\Documents\ArcGIS\Default.gdb\minus_ras2 with output format FGDBR.
On the other hand, if I explicitly designate raster objects using the format
distR = EucDistance(tmpFile, radius, 500) dist1 = radius - Raster("distR") dist2 = Raster("dist1") / radius
then I get the result
Runtime error Traceback (most recent call last): File "<string>", line 43, in <module> RuntimeError: ERROR 000732: Input Raster: Dataset distR does not exist or is not supported
Can anybody explain what the problem is and, most importantly, help me complete this code so it works properly?
Following up on Luke's suggestion, at this point I've now tried it with none of the files in a gdb, only the temporary files and the output raster in one, or all of the files in one (no spaces in the path/filenames). All files are in the same projection with the same extents. I've also tried it with 'Add results of geoprocessing operations to the display' both on and off.
I actually had it working the other night, but it crashed during operation and I must have changed something either in the program itself or in my ArcGIS configuration somewhere. This is more or less what the output is supposed to look like:
Now, however, whenever I enter the 'else' part of the 'if' statement (modified as below):
for pt in ptCursor: itemID = pt.getValue("RecordID") radius = pt.getValue("radius") query = "\"RecordID\" = " + str(itemID) arcpy.FeatureClassToFeatureClass_conversion(ptsLyr.name, os.path.split(tmpFile), os.path.split(tmpFile)[:-4], query) getmax = distR.maximum dist1 = Con(IsNull(distR),getmax,distR) dist2 = radius - dist1 dist3 = dist2 / radius if counter == 0: sumRas = dist3 else: sumRas = sumRas + dist3 sumRas.save(outfilename)
the sumRas raster gets replaced with no data, which obviously prevents me from accumulating a sum. I have no idea why this is happening and why I can't fix it. I've probably spent 40 hours on this now and I'm practically tearing my hair out.