I'm working on a little script that will do Local Density Analysis for
archaeology, and have run into an odd problem for a high-end GIS like GRASS.
I am having considerable difficulty in getting the count of points within a
given distance of another point. That is, for each point, I'd like to get
the number of points (in either the same or, better, another map layer)
within 100m, 200m, 300m, etc.
I first wrote my script using v.neighbors--reading the documentation rather
than trying it out (I know, bad idea)--followed by r.sum. It seemed perfect,
but it turns out that v.neighbors doesn't do what it seems it ought to do
given the documentation. When have time to figure out exactly what it IS
doing, I will see if I can clarify the docs.
R.neighbors is too restrictive in its window size (only 3-25) and the window
size is in grid units rather than in ground units (making my task
additionally complicated).
Creating a raster buffer with bands at the set of distances I need, followed
by r.statistics seemed like a good idea. But again, the actual behavior of
r.statistics does not seem to match quite what the docs imply--though the
docs also need to be updated a lot. Again, if I can be sure of what
r.statistics is doing I can try to update the docs.
Perhaps v.distance might be a place to start. However, I was trying to test
it on the Spearfish dataset to understand its behavior and stymied by the
lack of a dbf attribute file for any of the point data AND GRASS's apparent
inability to construct one from scratch (it seems like v.in.db will add
fields to an existing table but won't make one if it doesn't exist).
However, I'm still not sure what distances will end up being put into the
selected distance field if I choose the flag to calculate all distances
within a threshold. And I'm still left with a question about how to count
points within a set distance unless GRASS supports count(*) in its SQL
subset (it doesn't say that it does, but perhaps someone knows for sure).
As far as I can see, there are no distance functions in r.mapcalc.
Right now, I'm looking at doing a series of buffers, making each a mask and
running r.sum for each mask. This seems klunky but doable.
I am hoping that there is some kind of solution to this that I am missing.
Any suggestions? Thanks in advance.
Michael
____________________
C. Michael Barton, Professor of Anthropology
School of Human Evolution and Social Change
PO Box 872402
Arizona State University
Tempe, AZ 85287-2402
USA
Phone: 480-965-6262
Fax: 480-965-7671
www: <www.public.asu.edu/~cmbarton>