nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

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!

Vind meer berichten getagd met

Accepted answers

Anonymous User

I actually only need the ascii file as I can read this into R where I am doing my analysis.

WOW! That is good news...Well you can eliminate the last part then! The way I wrote that is to have it delete the ascii files after the conversion...I have eliminated the convert ascii to raster so you will just have ascii files. This will probably do all 45k points in a few short hours.

I have just attached a new toolbox, this one has 2 tools...The first one is the tool that converts each ascii to raster...You do not want to use this one. Use the (ascii) one. This will just write out ascii files and not convert to raster...These ascii files can still be used in ArcMap...If you pull them in they will be projected properly.

So if you use the (ascii) tool, it will hopefully be done before you leave the office today on all 45K points! I usually do not spend this much time working on a problem I find on the forums, but this one really interested me. I think this tool is useful. Let me know how it works for you!

One more note....Be sure to pull a few of the ascii rasters into ArcMap overlaid with your points to make sure they are centered...You also may want to use 110 meters for the distance parameter to be safe. Enjoy!

EucDistTools.zip

All comments

Anonymous User

No problem, can you send me say 20 points or so? I had it set up to do the OID field. I can change that as well...Use my personal email if you can share (caleb.mackey@gmail.com) and I can tailor it for your unique situation...This way we do not need to blow up this thread anymore, its getting pretty long!

And don't worry about wasting my time...This is a slow time at work for me. I think I will use this tool as well; this interests me.

KarenHornigold

Yes I have downloaded the new toolbox

.

I was already running the first script tool to create rasters for a subset of approx 600 points - it took 28 minutes!

But yes, they are shifted to the right. How can I correct this? I'm afraid to start messing about with your beautiful script!

Also the files came out with names 'dist_0', 'dist_1' etc., but I need them to take the VISITID code from the VISITID field in the point layer to be able to join this data with other datasets. This was in the very original script that you wrote. Would this be simple to do?

Sorry to take up more of your time!

Anonymous User

Yes...But first, be sure you download the new toolbox I just posted (probably as you were writting that question!).

So why I had to subtract 2 cell sizes is because it was using the original lower left corner because after my first run I saw that the Euclidean distance raster was off center by exactly 2 cell sizes...I guess I meant more by your coordinates, not by units (i.e. if coords are positive or negative). You may not have to change it to plus, I'm not sure. By subtracting that it centered them all for me. Just know that it shifts the cells to the correct position, lol. I hope that is a good enough answer.

KarenHornigold

One more note...The lines highlighted in blue may need to be changed to + instead of - depending on the units...I had to subtract 2 cell sizes to make sure it was centered.

Could you just explain this a bit more please? Why did 2 cell sizes need to be subtracted and what will affect this? My units will always be meters. I'd like to be sure I'm using the script correctly, and be able to check it has worked. Thanks!

Anonymous User

I actually only need the ascii file as I can read this into R where I am doing my analysis.

EucDistTools.zip

KarenHornigold

Wow, this is incredible! Thank you so much! I will cite your tool in any papers I submit.

I actually only need the ascii file as I can read this into R where I am doing my analysis. Currently I am reading in the rasters so I can just alter this line of the R script for the different file type! So do I just remove this part of the code:

   # ascii to raster
    Message('converting ascii files to raster...this may take a while')
    arcpy.env.workspace = ws
    for txt in arcpy.ListFiles('tmp*.txt'):
        number = int(''.join(n for n in txt if n.isdigit()))
        out = os.path.join(outws, 'dist_%s' %number)
        if arcpy.Exists(out): arcpy.Delete_management(out)
        arcpy.ASCIIToRaster_conversion(txt, out)
        arcpy.DefineProjection_management(out, SR)
        arcpy.Delete_management(txt)
        arcpy.Delete_management(txt[:-4] + '.prj')
        
    # delete templates and ascii folder
    for txt in arcpy.ListFiles('ascii_*'):
        arcpy.Delete_management(txt)
    os.rmdir(ws)

Anonymous User

Karen,

Scratch that, I have found a much improved and lightning fast way to do this (well lightning fast compared to the other way....)

Ok, this was driving me crazy so I spent several hours trying to figure out the most efficient way to do this...I think I finally found something that can do this as quickly as possible.

So, I tried to resolve this issue by creating a buffer around each point that is a little larger than your desired 100 unit neighborhood (I used 120). I then set the output extent environment to that of the buffer...This sped it up a little bit, but the euclidean distance still took a long time to run (since each raster was only 24x24 rows by columns). This got me thinking I could try making ascii files for each one. So I exported 5 or 6 of those rasters to ascii to look at each one...I found that every single ascii file was exactly the same except for the 3rd and 4th line (xllcorner, yllcorner). These coordinates are the x,y coordinates for the lower left corner for each raster. This matched the lower left corner for the extent of each buffer. This was good news...

So knowing that every ascii text file was exactly the same except for each point except for those 2 lines (because they all have the same distance neighborhood), I knew that I could take one of the ascii files and the accompanying prj file and just rename them and change the lower left corner coordinates.

So the script no just creates a buffer for each point and grabs the extent for that buffer. It then fixes the xllcorner/yllcorner coordinates (i.e. just moves a copy of that ascii grid to that lower left corner) and makes an ascii raster for each one....This part is very fast...The slow part will just be converting each ascii file to raster...but this is still about 50x faster than using the euclidean distance tool.

So how it works....It creates one euclidean distance raster for the first point as a template, then converts that raster to ascii and uses that ascii file to build all the others by just adjusting the lower left corner. Anyways it seems to work great...

I did 120 points in ~1 minute, and 1000 points in little over an hour...So this may work for all 45k of your points in a ~50 hours (seems to be a snowball effect on memory, still looking into that). If you have an i7 processor and 16GB of RAM you may have better luck than me; I have i5, 8GB. I'd say try it on 1000 points first and see how that goes. I don't know if you saw, but in post # 14 I added a script to divide the points by every 1000 (or however many you want). It may be better to run this in smaller chunks anyways to not kill your memory.

Here is the results window from running it as a script tool (with a lot of print statements removed of course)...The print statements can be taken out, they will only slow it down for a lot of points. I just left those in there from when I was testing.

Executing: DistanceFromEachPoint C:\TEMP\Points_divided\points_1000.shp C:\TEMP\Points_divided\rasters 100 FEET 10
Start Time: Fri Aug 09 09:36:16 2013
Running script DistanceFromEachPoint...
count: 1000 points...
buff_dist: 120.0
created euclidean distance template
Created ascii templates...
processed point ID: 0
processed point ID: 1
processed point ID: 2
processed point ID: 3
....
....
processed point ID: 998
processed point ID: 999
converting ascii files to raster...this may take a while
Complete (Elapsed time: 1:05:26.909)
Completed script DistanceFromEachPoint...
Succeeded at Fri Aug 09 10:41:43 2013 (Elapsed Time: 1 hours 5 minutes 26 seconds)

Here is the script, if you run by stand alone just change the parts highlighted in red...One note, I have this set up to work on the 'OID' field. You can change this to your 'VISITID' field. I just used OID as a default. Otherwise I have also attached a [bold] version 10.1 Script Tool[/bold]

'''
Euclidean distance for individual points-  This tool was written for cases where
  the user wants Euclidean distance rasters for each point in a feature class. This
  is useful when the neighborhoods will overlap.

Written by:  Caleb Mackey
Date:  8/8/2013
'''

import arcpy, os, shutil, sys
arcpy.CheckOutExtension('Spatial')
from arcpy.sa import *
from datetime import datetime as d
startTime = d.now()

try:
    odir = os.path.dirname(__file__)
except:
    odir = os.getcwd()

def Message(msg):
    print str(msg)
    arcpy.AddMessage(str(msg))
    
def bufferPointExtent(point_lyr, field, value, distance, units):
    query = '''{0} = {1}'''.format(field, value)
    arcpy.SelectLayerByAttribute_management(point_lyr, 'NEW_SELECTION', query)
    temp_buff = r'in_memory\buff'
    arcpy.Buffer_analysis(point_lyr, temp_buff, '{0} {1}'.format(distance, units))
    ext = arcpy.Describe(temp_buff).extent
    arcpy.Delete_management(temp_buff)
    return ext
    

def DistanceFromEachPoint(points, outws, distance, units, cellsize, *args):
    '''
                Params
    points:     points feature class
    outws:      output workspace for euclidean distance rasters
    distance:   distance for euclidean distance tool (units are same as points) i.e. 100
    units:      linear units (used for buffer) i.e. Feet or Meters

    '''
    # set output workspace
    arcpy.env.workspace = outws 
    arcpy.env.overwriteOutput = True
    arcpy.env.nodata = 'NONE'
    arcpy.env.pyramid = 'PYRAMIDS'
    SR = arcpy.Describe(points).spatialReference
    arcpy.env.outputCoordinateSystem = SR
    ws = os.path.join(odir,'ascii')
    if not os.path.exists(ws):
        os.makedirs(ws)

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

    # Make feature and count points
    lyr = arcpy.MakeFeatureLayer_management(points, 'points_lyr')
    count = int(arcpy.arcpy.GetCount_management(lyr).getOutput(0))
    Message('count: %s points...' %count)

    # generate template ascii file/prj file
    buff_dist = float(distance) + (float(cellsize)*2)
    buff_ext = bufferPointExtent(lyr, oid, min(ids), buff_dist, units)
    arcpy.env.extent = buff_ext
    query = '''{0} = {1}'''.format(oid, min(ids))
    arcpy.SelectLayerByAttribute_management(lyr, 'NEW_SELECTION', query)
    sample = EucDistance(lyr, distance, cellsize)
    eucDist = os.path.join(odir, 'eucDist_tmplt')
    if arcpy.Exists(eucDist): arcpy.Delete_management(eucDist)
    sample.save(eucDist)
    Message('created euclidean distance template')
    arcpy.env.extent = None # if this is NOT set to none, will cause no selected features

    # templates
    ascii = os.path.join(ws, 'ascii_template.txt')
    arcpy.RasterToASCII_conversion(eucDist, ascii)
    arcpy.Delete_management(eucDist)
    prj = ascii[:-4] + '.prj'
    Message('Created ascii templates...')

    # Get ascii lines
    lines = open(ascii,'r').readlines()

    # Create ascii raster for each point 
    for i in ids:
        ext = bufferPointExtent(lyr, oid, i, distance, units)
        if ext: 
            lines[2] = 'xllcorner     %s\n' %str(float(str(ext.lowerLeft).split()[0]) - (float(cellsize)*2))
            lines[3] = 'yllcorner     %s\n' %str(float(str(ext.lowerLeft).split()[1]) - (float(cellsize)*2))
            copy = os.path.join(ws, 'tmp_%s.txt' %i)
            prj_copy = shutil.copy(prj, copy[:-4] + '.prj')
            with open(copy, 'w') as w:
                w.writelines(lines)
            Message('processed point ID: %s' %i)

    # ascii to raster
    Message('converting ascii files to raster...this may take a while')
    arcpy.env.workspace = ws
    for txt in arcpy.ListFiles('tmp*.txt'):
        number = int(''.join(n for n in txt if n.isdigit()))
        out = os.path.join(outws, 'dist_%s' %number)
        if arcpy.Exists(out): arcpy.Delete_management(out)
        arcpy.ASCIIToRaster_conversion(txt, out)
        arcpy.DefineProjection_management(out, SR)
        arcpy.Delete_management(txt)
        arcpy.Delete_management(txt[:-4] + '.prj')
        
    # delete templates and ascii folder
    for txt in arcpy.ListFiles('ascii_*'):
        arcpy.Delete_management(txt)
    os.rmdir(ws)
    Message('Complete (Elapsed time: %s)' %(str(d.now() - startTime)[:-3]))

if __name__ == '__main__':

    try:
        # always try to run as script tool first, otherwise, exception will allow stand alone params
        # Get params
        argv = tuple(arcpy.GetParameterAsText(i) for i in range(arcpy.GetArgumentCount()))
        DistanceFromEachPoint(*argv)

    except:
        arcpy.AddMessage('Script Tool failed')
        arcpy.AddError(arcpy.GetMessages(2))
        # stand alone
        points = r'C:\TEMP\Stream_points.shp'
        work = r'C:\TEMP\Output_Rasters\ascii_test\Rasters'

        # run it
        DistanceFromEachPoint(points, work, 100, 'Feet',  10)

One other option...If you don't mind them being in ascii format, it would be done in less than 2 hours probably. I am not sure if ascii rasters are valid inputs for processing tools though. The conversion of these to raster is what slows the tool down; but is still faster than running euclidean distance for each one.

One more note...The lines highlighted in blue may need to be changed to + instead of - depending on the units...I had to subtract 2 cell sizes to make sure it was centered. Anyways, this worked perfectly for me, and much faster! Enjoy!

EucDistTools.zip

KarenHornigold

Wow, this is amazing thank you! You're right, I have unnecessary no data cells for my whole extent. Thanks for your time and effort, I don't know what I would have done without your help!

Interestingly, yesterday I reran the original script on data from one county (560 points) and it was racing through the script- it took only 30 minutes! The things I changed were - 1) output raster format to .img and 2) scratch workspace from default geodatabase(! - don't know why that was the default) to c:\temp. It was really quick even though the extent was still for the whole country.

Nevertheless I would still like to use your buffer version to reduce file size, but I got the following error:

Runtime error 
Traceback (most recent call last):
  File "<string>", line 32, in <module>
AttributeError: DescribeData: Method extent does not exist

Thanks again

Anonymous User

I was rereading this post, and I saw that after you ran the first code I posted on 10 points, you said it took ten minutes...I thought that was extremely slow for how small each raster is. So I did a quick test on 36 points and it took way too long as well...Then I saw the size of each raster (600 MB). So from there I knew something was definitely wrong since these rasters are "supposed" to be small. I figured out that each raster is actually very large due to each raster filling in "NODATA" values for the entire extent of the points feature class.....so in your case, the extent is probably massive so there are a bunch of extra "NODATA" cells being created.

So a workaround:

I added a step to the original code to create a buffer (120 m) around each selected point...I then grab the extent of that buffer and set the "extent" environment to the extent of that buffer so now instead of being extremely large rasters they are now 24x24 and only 2.25 KB. This works much faster now. It will still take a very long time but it now does all 36 of my points in ~3 as opposed to ~ 12 minutes. The best performance gain here is it will take up a lot less space on disk...So try this on say 100 points and see how long it takes. If it it does not take too long, I would say try it on your entire data set and *hopefully* it will finish in a couple days.

import arcpy
from datetime import datetime as d
startTime = d.now()
arcpy.CheckOutExtension('Spatial')
from arcpy.sa import *

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

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

# clean up
for raster in arcpy.ListRasters('dist_*'):
    arcpy.Delete_management(raster)

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

# 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)
    temp_buff = 'buff_%s.shp' %i
    arcpy.Buffer_analysis(lyr, temp_buff, '120 METERS')
    extent = arcpy.Describe(temp_buff).extent
    arcpy.env.extent = extent
    eucDist = EucDistance(lyr, 100, 10)
    eucDist.save('dist_%s' %i)
    arcpy.Delete_management(temp_buff)
    print 'Created Euclidean Distance raster for point ID: %s' %i
print 'Done'
print '(Elapsed time: %s)' %(str(d.now() - startTime)[:-3])

You can leave out the "clean up" part if you want. This is still pretty slow, but should be at least 30x faster than before.

I did just spot some interesting behavior...About 1/4 of the rasters are still coming out with the large extent (of the whole point fc)....There could be something going on behind the scenes here or a bug...They should ALL come out the extent of the buffer...I think I will contact tech support and see if they can figure out the issue here. In the meantime, you could try the above code on a decent amount of points (i.e. 100) and see how long it takes. I added a statement to test the time.

Anonymous User

More important than this is the processing time issue. Is there a way I could batch process subsets of my points file and run them simultanously outside of ArcGIS with a Python programme?

I'm not sure if this will be any faster, but since you have ArcGIS 10.1 you can download Service Pack 1 (if you have not done so already) and then download the 64" rel="nofollow" target="_blank">http://resources.arcgis.com/en/help/main/10.1/index.html#/Background_Geoprocessing_64_bit/002100000040000000/]64 bit background processing. This will allow you to utilize the 64 bit python as a stand alone script.

The 64 bit processing does not necessarily increase speed, but it does take advantage of more processing power.

As for why your parameters are getting switched around, maybe you did not set up the script tool correctly? Based on your post reporting the errors, it looks like it the params were set up correctly but I am not sure as I cannot see your actual tool. The script looks fine to me though (after the cursor change).

I would recommend downloading the 64 bit background processing and run it as a stand alone script. ( The 64 bit BP can be found in your customer support page under downloads) This will take a very long time to run either way.

One last idea,

how should I change this line

Code:

query = '{0} = {1}'.format(arcpy.AddFieldDelimiters(points,'VISITID'),i)

to change VISITID to the field selected via the tool dialogue box?

You had it correct in your last post to be interpreted from the tool dialog (get param as text). So if you are using the "field" variable, just leave it as you had it:

query = '{0} = {1}'.format(arcpy.AddFieldDelimiters(points,field),i)

as a stand alone script you can just use 'VISITID' or make a field variable for it. Your choice.

By the way...You mentioned you wanted to split the data set up into smaller chunks? You could use this script to split it into groups of 1000 (or however many) and run the first code on one or 2 sets of data a day. This shows how you can split your points into 1000's:

import arcpy, itertools arcpy.env.workspace = r'C:\TEMP\Points_divided' arcpy.env.overwriteOutput = True points = r'C:\TEMP\Address_points.shp' oid = arcpy.AddFieldDelimiters(points, arcpy.Describe(points).oidFieldName) ids = list(r[0] for r in arcpy.da.SearchCursor(points,['OID@'])) sql_dict = {} num = 1000 #splits by every 1000 mini = 0 for key, group in itertools.groupby(ids, lambda k: k//num): maxi = int(max(group)) sql_dict[key] = [mini, maxi] mini += num # Make feature layer lyr = arcpy.MakeFeatureLayer_management(points, 'points_lyr') for v in sql_dict.values(): query = '{0} >= {1} AND {2} <= {3}'.format(oid, v[0], oid, v[1]) print query arcpy.SelectLayerByAttribute_management(lyr, 'NEW_SELECTION', query) arcpy.CopyFeatures_management(lyr, 'points_%s.shp' %int(v[1] +1)) print 'Created points_%s.shp' %int(v[1] + 1) print 'done'

KarenHornigold

One last idea,

how should I change this line

query = '{0} = {1}'.format(arcpy.AddFieldDelimiters(points,'VISITID'),i)

to change VISITID to the field selected via the tool dialogue box?

KarenHornigold

Now it creates an output but the raster values are strange (they are not the same as when I run the script with hard-coded paths) and it also doesn't save to the location I am specifying, so I'm still doing something wrong! I may give up on creating a script tool with this.

More important than this is the processing time issue. Is there a way I could batch process subsets of my points file and run them simultanously outside of ArcGIS with a Python programme?

Anonymous User

Well that is interesting....I don't know why you are getting that error. That means it does not know the format of the output raster. I deliberately did not add a file extension to the name of the raster so it would be the Esri GRID format. I usually just use GRID or tif rasters.

Maybe try and save at as a .tif?

eucDist.save('dist_%s.tif' %i)

Curtis Price is the raster master...perhaps he can weigh in on this?? He may even have a better methodology for this whole process.

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)

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)

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#//018w00000011000000]da 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

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

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.

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

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.

KarenHornigold

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

With gratitude

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.

Doe mee aan het gesprek

Aangemelde leden kunnen berichten plaatsen, updates volgen en meer. Nieuw hier? Registreer een gratis account.

Snelle Links

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.