Iterating through point layer to create & sum distance rasters

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 Filesproject.mxd"
ptsFile = r"C:Project Filesdatapoints.shp"
sumFile = r"C:Project FilesdataProject.gdbsumbuffer"
tmpFile = r"C:Project Filesdataone_point.shp"

mxd = amap.MapDocument(mapFile)
df = amap.ListDataFrames(mxd)[0]
mlayers = amap.ListLayers(mxd)
for lyr in mlayers:
    if lyr.name == os.path.split(ptsFile)[1][:-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)[0], os.path.split(tmpFile)[1][:-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:UsersJohnRDocumentsArcGISDefault.gdbminus_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)[0], os.path.split(tmpFile)[1][:-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.