Skip to content
Prev 3071 / 29559 Next

sp basics !

Hi Didier,
On Mon, 28 Jan 2008, Didier Leibovici wrote:

            
The most recent public ones are from the Imperial College course:

http://www.bias-project.org.uk/ASDARcourse/
(where this isn't geoR the package? - to avoid confusion, maybe 
R-spatial?)
Is this distance as in Euclidean, or is it Great Circle? How precise do 
you need to be? Will you accept the boundary points of one as points to 
the boundary of the other as a line? Have a look at the spatstat internal 
function distppll(), which is used there for point pattern analysis within 
a window, but works well for this version of points to line and Euclidean 
distance. Watch the l= argument, it is in segments and has x0, y0, x1, y1 
columns for each segment. This is messy, but:

library(spdep)
example(eire)
# to load a small data set, Irish counties

eire$all <- rep(1, 26)
eire_per <- unionSpatialPolygons(eire, eire$all)
# to make a perimeter to measure to

sapply(slot(eire_per, "polygons"), function(x) length(slot(x,
   "Polygons")))
sapply(slot(eire_per, "polygons"), function(x) sapply(slot(x, "Polygons"),
   slot, "area"))
# to see which bit of Co. Mayo is the island - since it is the second,
# choose [[1]]

ll <- slot(slot(slot(eire_per, "polygons")[[1]], "Polygons")[[1]],
   "coords")
# get the perimeter coordinates

pls <- slot(eire, "polygons")
lgths <- sapply(pls, function(x) length(slot(x, "Polygons")))
sapply(slot(pls[[which(lgths == 2)]], "Polygons"), slot, "area")
# to see which bit of Co. Mayo is the island, same diagnosis, we
# can use [[1]] in the work loop

lll <- cbind(ll[-nrow(ll),,drop=FALSE], ll[-1,,drop=FALSE])
# to set up the segments from the perimeter coordinates

res <- vector(mode="list", length=length(pls))
for (i in seq(along=res)) {
   crds <- slot(slot(pls[[i]], "Polygons")[[1]], "coords")
   res[[i]] <- apply(distppll(p=crds, l=lll), 1, min)
}
# actually do the work
res

gets you the shortest distances to each point on the boundary of each 
county. Note that Donegal has 0 distance for all points, because there are 
no boundary points along its land border, though there could have been.
Using distppll() here is low level, so it just loops over points for each 
segment (or vice-versa), not picky about topology. But the S4 
SpatialPolygons are easier to chuck around, because we know definitely 
what is inside them. "polylist" objects were S3.
No, though if someone felt like trying out an OGR virtual driver and GEOM, 
there are robust methods for classes out there. Or you could push the sata 
out to PostGIS, modern MySQL, or any supported database and Terralib using 
the aRT interface. At least one of these may support geographical 
coordinates, i.e. Great Circle distances.

If you make any progress with alternatives, please let us know.

Hope this helps,

Roger