[R-sig-Geo] sp basics !

[Roger Bivand]

Hi Didier,

On Mon, 28 Jan 2008, Didier Leibovici wrote:

> Hi all,
> I am still novice in using sp and/or Maptools and  I find difficult to
> locate basic tutorials and general description ...
> I am probably some good source of info?
>

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

http://www.bias-project.org.uk/ASDARcourse/

> My problem is very simple and I brlieve it covers some basic handling of
> "GIS" data desirable in geoR:

(where this isn't geoR the package? - to avoid confusion, maybe 
R-spatial?)

>
> I would like to read some shapefiles of polygons( I know how to do that)
> and loop through all the polygons to be able to compute the shortest
> distance to one another (fix) polygon
> e.g. distance to the sea or to a montain

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.

> -I had at one point a polylist but probably I need to have "consistant"
> polygon for each element of the list ?

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.

> -is there any distance method supported by a class "polygon" which can
> compute the distance to another geometry?

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

>
> thanks
>
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no