Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have a point file, hourly relocations of an animal, and I want to be able to place a buffer around each point and calculate the number of subsequent points which are within the buffer zone. I am looking for a method which will work along the point file, like a moving window, which will only count those subsequent points which are within the buffer zone.

For instance at point 10 I place a buffer of 1500m and I want to know whether point 11 is within the buffer and if so then whether point 12 is within the buffer and so forth. I don't want to know whether point 52 is within the buffer zone unless all of the previous points have been within the buffer. I also don't want to know whether points 9 or 8 etc. are within the buffer.

I have found and tried an additional toolbox called "moving window point analysis" which works as a moving window on the point file. This works well, but very slowly, and includes all points that are within the buffer zone even if they are not consecutive points. I can't find a way to just make it look at consecutive points.

I would like a method which will provide an output table as I have a lot of data points to look at in this way.

I am using ArcGIS 10. Any help that anyone can provide would be greatly appreciated.

Thanks

share|improve this question
    
Your points likely originated as a list of (x,y,time) data. Would you be open to pre-processing these data (outside ArcGIS) to obtain the desired information? –  whuber Feb 17 '12 at 16:03
    
If that makes it easier then definitely. I am also processing the data using AdehabitatLT in R which calculates distance and bearings etc. I understand the process suggested by Sylvester below but I am struggling to know where to start as I'm not really sure which tools I need to use. –  James Feb 17 '12 at 16:28
    
Ah! Since you're already using R, let's explore R-based solutions then. –  whuber Feb 17 '12 at 17:07
    
There is a sliding window function within AdehabitatLT "sliwinltr" but I don't know how to use it in this instance. I don't even know if it can be used in this way. –  James Feb 17 '12 at 17:12
add comment

3 Answers

Given a list of point locations (preferably in projected coordinates, so that distances are easy to compute), this problem can be solved with five simpler operations:

  1. Compute point-point distances.

  2. For each point i, i = 1, 2, ..., identify the indexes of those points at distances less than the buffer radius (such as 1500).

  3. Restrict those indexes to be i or greater.

  4. Retain only the first consecutive group of indexes having no break.

  5. Output the count of that group.

In R, each of these corresponds to one operation. To apply this sequence to each point, it's convenient to encapsulate most of the work within a function we define, thus:

#
# forward(j, xy, r) counts how many contiguous rows in array xy, starting at index j,
#                   are within (Euclidean) distance r of the jth row of xy.
#
forward <- function(j, xy, r) {
  # Steps 1 and 2: compute an array of indexes of points within distance r of point j.
  i <- which(apply(xy, 1, function(x){sum((x-xy[j,])^2) <= r^2}))
  # Step 3: select only the indexes at or after j.
  i <- i[i >= j]
  # Steps 4 and 5: retain only the first consecutive group and count it.
  length(which(i <= (1:length(i) + j)))
}

(See below for a more efficient version of this function.)

I have made this function flexible enough to accept various point lists (xy) and buffer distances (r) as parameters.

Normally, you would read a file of point locations (and, if necessary, sort them by time). Here, to show this in action, we will just generate some sample data randomly:

# Create sample data
n<-16                                     # Number of points
set.seed(17)                              # For reproducibility
xy <- matrix(rnorm(2*n) + 1:n, n, 2) * 300
#
# Display the track.
plot(xy, xlab="x", ylab="y")
lines(xy, col="Gray")

Figure

Their typical spacing is 300*Sqrt(2) = about 500. We do the calculation by applying this function to the points in the array xy (and then tacking its results back on to xy, because this would be a convenient format for export to a GIS):

radius <- 1500
z <- sapply(1:n, function(u){forward(u,xy,radius)})
result <- cbind(xy, z)                              # List of points, counts

You would then further analyze the result array, either in R or by writing it to a file and importing it into other software. Here is the result for the sample data:

                        z
  [1,]   -4.502615  551.5413 4
  [2,]  576.108979  647.8110 3
  [3,]  830.103893 1087.7863 4
  [4,]  954.819620 1390.0754 3
...
 [15,] 4977.361529 4146.7291 2
 [16,] 4783.446283 4511.9500 1

(Remember that the counts include the points at which they are based, so that each count must be 1 or greater.)


If you have many thousands of points, this method is too inefficient: it computes far too many point-to-point distances that are unnecessary. But because we have encapsulated the work within the forward function, the inefficiency is straightforward to fix. Here is a version that will work better when more than a few hundred points are involved:

forward <- function(j, xy, r) {
   n <- dim(xy)[1]     # Limit the search to the number of points in xy
   r2 <- r^2           # Pre-compute the squared distance threshold
   z <- xy[j,]         # Pre-fetch the base point coordinates
   i <- j+1            # Initialize an index into xy (just past point j)

   # Advance i while point i remains within distance r of point j.
   while(i <= n && sum((xy[i,]-z)^2) <= r2) i <- i+1

   # Return the count (including point j).
   i-j
}

To test this, I created random points as previously but varied two parameters: n (the number of points) and their standard deviation (hard-coded as 300 above). The standard deviation determines the average number of points within each buffer ("average" in the table below): the more there are, the longer this algorithm takes to run. (With more sophisticated algorithms the run time won't depend as much on how many points are in each buffer.) Here are some timings:

 Time (sec)    n    SD  Average  Distances checked per minute
 1.30       10^3     3  291      13.4 million
 1.72       10^4    30   35.7    12.5
 2.50       10^5   300    3.79    9.1
16.4        10^6  3000    1.04    3.8
share|improve this answer
    
This looks like the perfect solution. However the code doesn't run from "z <- sapply(1:n), function(u){forward(u,xy,radius)})" as it says: "unexpected ',' in "z <- sapply(1:n)," If I remove the comma it then says: Error: unexpected 'function' in "z <- sapply(1:n) function" Any ideas why this might be so? –  James Feb 17 '12 at 18:21
    
Sorry; there's a typo there: I'll remove the extraneous ")". (I made a last minute simplification of this code. It's been tested more times than I care to admit!) –  whuber Feb 17 '12 at 18:25
1  
Thats great, it's running now. I just loaded my data in as xy for simplicity to check it works. It takes a little time to run as you mentioned but seems to have done it all correctly. I will manually double check a few with my GIS map but it looks good so far. Thanks for helping me work this out, very keen to learn in both GIS and R and i'm on the steep learning curve. –  James Feb 17 '12 at 18:42
1  
I edited the reply to provide a solution with greatly improved scalability. It is now capable of handling paths containing millions of points. –  whuber Feb 19 '12 at 16:25
1  
I ran the original code with point files of 2000 entries which took a couple of hours each, as you say there were far too many unused point to point locations extending the processing time. The edit above looks like a neat solution and I will try this on the same data and see how much faster it is. Thanks for the effort with producing and editing the function. –  James Feb 20 '12 at 9:58
show 3 more comments

I think your best bet is to script a little routine using ArcPy. I'd create something like this pseudo-code:

Select all points sorted by location-id    
Iterate for each point:
    Select points by location using a distance (no need to create buffer geometry)
    Sort points by location-id
    Set a variable to the value of your reference id
    consecutive-counter = 0
        Iterate your selection:
            Is the location-id of the first (or next) point equal to variable + 1?
            if 'yes' increment consecutive-counter by 1
            Repeat until 'no'

I'm not sure what you want to do with the information but I suppose you could create a field in your table and update it with the count of consecutive points (if so, add the field first).

I would recommend making feature layers (like a database table view but for features in Arc). Make two from the original data and open an update cursor on the first set specifying your over-all sort (because ESRI doesn't honour full SQL queries). Use the second one to select by location from and open a Search Cursor on the resultant selection-set.

[EDIT AS PER Jame's REQUEST] Rough it out using Model Builder. If you have not used Model Builder before, all you have to do is right-click in arcToolbox. Select 'Add Tool box'. The right-click on the new toolbox and click 'New->model'. Once you have a new model window, drag and drop the tools and data you need into the window and visually link them together (using the little arrow tool). When you have got as far as you can (you'll not be able to add your cursors here), use the option in the Model Builder's File Menu to export to Python. That will get you most of the way there. It is auto-generated code so will be a bit nasty but functional. Then use the links in my answer above to understand and add the cursors.

If you are new to Python, don't be frightened of writing code! Python is a very easy scripting language to get results from quickly. Esri has some guidance on it too.

If you get stuck with your code, post it to this forum and ask for help. There's a LOT of people here that can help. One warning - make sure you use the right help from ESRI. They massively changed the Python API between versions 9.x and 10 (respectively arcgisscripting and arcpy). So if you are using ArcGIS 9.x find the equivalent links to mine!

share|improve this answer
    
This looks like what I would like to do. However, I am not currently using code within ArcGIS, merely selecting from pre-defined options. How would I begin to use/generate the suggested method above? I would like the output to be either a new table or a new field added to the table with the count of consecutive points. –  James Feb 17 '12 at 15:28
    
See my edit to my main post. –  MappaGnosis Feb 17 '12 at 15:43
    
JTB, please log on using the same account you used to post this question so that you can post comments. (To make it easier, I merged the JTB account with the James account.) –  whuber Feb 17 '12 at 15:50
    
Sorry about the account change. I posted the original question as a new user but then I couldn't get back into that account as I didn't have a password etc. Therefore I created another account JTB which i will use from now on (hopefully). I will start the model builder as suggested by Sylvester, but having never used it before it may take me a little time to get used to it and to work out what tools etc to use. I will come back with progress and questions. Thanks –  James Feb 17 '12 at 16:07
    
Sylvester - I think I understand the process but I am at a loss to know which tool(s) I really need to start with. Distance? Buffer? Near? I don't even know if there is a correct tool for this problem which will do what has been mentioned above. I am keen to learn but am very much at the beginning. –  James Feb 17 '12 at 16:32
add comment

You can use model builder in ArcGIS to find consecutive ID values. I exported my model as a python script. The code will generate a new shp that has consecutive ID values. !ID! is the base ID field. You will have to update the point2.shp path, name, and ID field name to match your case.

# Import arcpy module
import arcpy


# Local variables:
point2_shp = "C:\\temp\\point2.shp"
point2_shp__2_ = "C:\\temp\\point2.shp"
point2_shp__4_ = "C:\\temp\\point2.shp"
Freq_dbf__2_ = "C:\\temp\\Freq.dbf"
point2_shp__5_ = "C:\\temp\\point2.shp"
point2__2_ = "C:\\temp\\point2.shp"
point2__4_ = "C:\\temp\\point2.shp"
Freq_dbf = "C:\\temp\\Freq.dbf"
PointConsecutive_shp = "C:\\temp\\PointConsecutive.shp"

# Process: Add Field
arcpy.AddField_management(point2_shp, "AUTOID", "LONG", "", "", "", "", "NON_NULLABLE", "NON_REQUIRED", "")

# Process: Calculate Field
arcpy.CalculateField_management(point2__2_, "AUTOID", "autoIncrement()", "PYTHON_9.3", "rec=0\\ndef autoIncrement():\\n global rec\\n pStart = 1 #adjust start value, if req'd \\n pInterval = 1 #adjust interval value, if req'd\\n if (rec == 0): \\n  rec = pStart \\n else: \\n  rec = rec + pInterval \\n return rec\\n")

# Process: Add Field (2)
arcpy.AddField_management(point2_shp, "DIFF", "LONG", "", "", "", "", "NON_NULLABLE", "NON_REQUIRED", "")

# Process: Calculate Field (2)
arcpy.CalculateField_management(point2__4_, "DIFF", "!ID! - !AUTOID!", "PYTHON_9.3", "")

# Process: Frequency
arcpy.Frequency_analysis(point2_shp__2_, Freq_dbf, "DIFF", "")

# Process: Join Field
arcpy.JoinField_management(point2_shp__4_, "DIFF", Freq_dbf__2_, "DIFF", "")

# Process: Select
arcpy.Select_analysis(point2_shp__5_, PointConsecutive_shp, "\"FREQUENCY\" >1")
share|improve this answer
    
I am not sure what the code above does and how it answers my query. I appreciate the help but this is all quite new to me having never used Python or model builder. Do I change the ID field for each process listed to the ID in the dataset? –  James Feb 17 '12 at 17:36
    
@James, It depends if your ID changes. To use the code just copy and paste the code above and save it to a blank txt file. Update the code to match your point.shp name and path. Then, change the ID name in the Calculate Field (2) code section to match your point.shp ID field. Save the txt file, and in windows explorer re-name the file with a .py extension. Right click the file and open with python.exe to test. –  artwork21 Feb 17 '12 at 18:54
    
Ideally, this script could be plugged into a script tool, which also handles the buffering and feature selection. help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/… –  artwork21 Feb 17 '12 at 18:59
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.