Loop through features and create euclidean distance raster per feature

[KarenHornigold]

Hi,

I am completely new to Python, but need to turn to it due to the nature of my research.

I have a point shapefile with 45,000 points ('visits') spread across the UK and I want to obtain a euclidean distance raster (with a cut off of say 100km).

I know the tool I need to use to create the distance raster is Euclidean Distance, but I need to iteratively select each point (by 'visitID') from the shapefile and create a distance raster for each one in turn.

Can anyone direct me to useful resources or share their scripting wisdom?

Thank you in advance!

[Anonymous User]

Something like this would work, but it would be painfully slow to make a separate distance raster for all 45k points. This selects each point one at a time by its 'visitID' and from the selection runs a Euclidean Distance with max dist set for 100. I used a cell size of 10. Each raster is named "dist_" + the visitID of the point.

import arcpy
arcpy.CheckOutExtension('Spatial')
from arcpy.sa import *

# set output workspace
arcpy.env.workspace = r'C:\TEMP\Output_Rasters'
arcpy.env.overwriteOutput = True

# your points
points = r'C:\TEMP\Stream_points.shp'

# oid field and list of id's 
ids = list(r.visitID for r in arcpy.SearchCursor(points))

# loop through table to make a distance raster for each point
# EucDIstance(points, max_distance, cell_size)  <-- parameters
lyr = arcpy.MakeFeatureLayer_management(points, 'points_lyr')
for i in ids:
    query = '{0} = {1}'.format(arcpy.AddFieldDelimiters(points,'visitID'),i)
    arcpy.SelectLayerByAttribute_management(lyr, 'NEW_SELECTION', query)
    eucDist = EucDistance(lyr, 100, 10)
    eucDist.save('dist_%s' %i)
    print 'Created Euclidean Distance raster for point ID: %s' %i
print 'Done'

Make sure the 'visitID's highlighted in red are the correct field name and matching case.

[KarenHornigold]

Thank you for your quick reply and for providing code! I shall try it out and let you know.

With gratitude

[KarenHornigold]

So I just ran the script on a sample of 10 points and it works perfectly.

However, it took 10 minutes for 10 points, so a quick scaling up led to an estimate of one month to run on 45,000 points!

Is there any alternative, more efficient way?

Thanks again.

[Anonymous User]

Is there any alternative, more efficient way?

Sure, running Euclidean Distance on the whole shapefile at once!  Of course, you do not always have this luxury.  Is it absolutely necessary to have an individual raster for each point?  As long as your points aren't within 100 m of each other, there should not be any overlap. 

If you don't care so much about having ugly rasters you could use a cell size of 20.  That could speed it up a little bit, but it will still take forever with 45k points.

[KarenHornigold]

Is it absolutely necessary to have an individual raster for each point?

Yes, I need one raster per point as I want to find out distance weighted population surrounding each point. Also my points are very dense in some places with only a few meters between, so there would be an overlap whatever distance cut off I decide to use.

[KarenHornigold]

I attempted to turn the script into a custom tool, but can't get it working. The problem seems to be with how to replace 'VISITID' as a specified field, to use instead the field that is selected in the dialogue box. Can you take a look at my script to see what I'm doing wrong?

import arcpy
from arcpy import env
arcpy.CheckOutExtension('Spatial')
from arcpy.sa import *

# set parameters (needed for dialgue box)
points = arcpy.GetParameterAsText(0)
field = arcpy.GetParameterAsText(1)

# set output workspace
outWorkspace = arcpy.GetParameterAsText(2)
arcpy.env.overwriteOutput = True

# oid field and list of id's
ids = list(r.field for r in arcpy.SearchCursor(points))

# loop through table to make a distance raster for each point
# EucDIstance(points, max_distance, cell_size)  <-- parameters
lyr = arcpy.MakeFeatureLayer_management(points, 'points_lyr')
for i in ids:
    query = '{0} = {1}'.format(arcpy.AddFieldDelimiters(points, field),i)
    arcpy.SelectLayerByAttribute_management(lyr, 'NEW_SELECTION', query)
    eucDist = EucDistance(lyr, 10000, 1000)
    eucDist.save('dist_%s' %i)
    print 'Created Euclidean Distance raster for point ID: %s' %i
print 'Done'

Thanks

[Anonymous User]

Since you are wanting to supply the field name as a variable, try changing this line:

# oid field and list of id's 
ids = list(r.getValue(field) for r in arcpy.SearchCursor(points))

Note that the above is the 10.0 cursor. I was not sure if you have 10.1 or 10.2 yet, but if you do, you can use a da" rel="nofollow" target="_blank">http://resources.arcgis.com/en/help/main/10.1/index.html#//018w00... cursor (data access) like this:

# oid field and list of id's 
ids = list(r[0] for r in arcpy.da.SearchCursor(points, [field]))

So if you have 10.0 use the first option, or if you have 10.1 or above use the second. Also, for future reference, try posting your code using the code tags by using the # button or you can manually type it like this:

[noparse]

if condition:
     # do stuff
# more code

[/noparse]

By using the code tags, it will make the code look like this:

if condition:
     # do stuff
# more code

Indentation is very important in python so it makes it easier to spot errors when you use the code tags. You can see the full BB Code List [url=http://forums.arcgis.com/misc.php?do=bbcode]here[/url] Good luck!

[KarenHornigold]

Thanks for your reply!

I got the following error message:

Messages
Executing: IterativeEuclideanDistance C:\Docs\ArcGIS\Dorset_sample.shp VISITID C:\Docs\ArcGIS\DORSET
Start Time: Thu Aug 08 14:25:23 2013
Running script IterativeEuclideanDistance...

Traceback (most recent call last):
  File "C:\Docs\ArcGIS\Scripts\IterativeEucDist.py", line 15, in <module>
    ids = list(r[0] for r in arcpy.SearchCursor(points,[field]))
  File "C:\Docs\ArcGIS\Scripts\IterativeEucDist.py", line 15, in <genexpr>
    ids = list(r[0] for r in arcpy.SearchCursor(points,[field]))
TypeError: 'Row' object does not support indexing

Failed to execute (IterativeEuclideanDistance).
Failed at Thu Aug 08 14:25:23 2013 (Elapsed Time: 0.00 seconds)

[KarenHornigold]

Sorry I noticed I missed the 'da' in searchcursor. So I ran it again with this corrected and got the following error:

Messages
Executing: IterativeEuclideanDistance C:\Docs\ArcGIS\Dorset_sample.shp VISITID C:\Docs\ArcGIS\DORSET
Start Time: Thu Aug 08 14:32:28 2013
Running script IterativeEuclideanDistance...

Traceback (most recent call last):
  File "C:\Docs\ArcGIS\Scripts\IterativeEucDist.py", line 24, in <module>
    eucDist.save('dist_%s' %i)
RuntimeError: ERROR 010093: Output raster format UNKNOWN is unsupported.

Failed to execute (IterativeEuclideanDistance).
Failed at Thu Aug 08 14:33:10 2013 (Elapsed Time: 42.00 seconds)