# Iterating through a 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 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)[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")
# 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:\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)
``````

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?

2 April

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")
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)

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.

You could try setting the output workspace with `arcpy.env.workspace` to be a folder not a file GDB.