Message-ID: <alpine.LRH.2.00.1102162140140.29370@reclus.nhh.no>
Date: 2011-02-16T20:56:56Z
From: Roger Bivand
Subject: classify point by nearest polygon
In-Reply-To: <1297816552305-6030040.post@n2.nabble.com>
On Tue, 15 Feb 2011, rundel wrote:
>
> This can also be done with rgeos, a simple example is below. The 2nd half is
> relevant to what you are trying to accomplish the rest is just generating
> lakelike polygons.
>
Maybe it was just my mailer breaking lines at the wrong places, but a
couple of updates are needed for this to fly:
>
set.seed(1)
txt = rep("",5)
for(i in 1:5) {
xc = floor(runif(1,1,25))
yc = floor(runif(1,1,25))
xpts = round(rnorm(10,0,0.5),3)
ypts = round(rnorm(10,0,0.5),3)
# txt[i] = paste("MULTIPOINT(",paste( "(",xpts+xc,"
# ",ypts+yc,")",sep="",collapse="," ),")",sep="")
tmp <- paste("(", xpts+xc, " " ,ypts+yc,")", sep="", collapse=",")
txt[i] = paste("MULTIPOINT(", tmp, ")", sep="")
}
wkt = paste( "GEOMETRYCOLLECTION(",paste(txt,collapse=","),")", sep="")
z2=gConvexHull(readWKT(wkt),byid=TRUE)
lakes = gBuffer(z2,byid=TRUE,width=1,quadsegs=10)
plot(lakes,col=1:5+1)
pts = matrix( runif(100,1,25),ncol=2)
#pts_wkt = paste( "GEOMETRYCOLLECTION(",
#paste(paste("POINT(",pts[,1],"
#",pts[,2],")",sep=""),collapse=",") ,")", sep="")
tmp <- paste("POINT(",pts[,1]," ",pts[,2],")",sep="")
pts_wkt = paste("GEOMETRYCOLLECTION(", paste(tmp, collapse=",") ,")",
sep="")
pts_sp = readWKT(pts_wkt)
#cols = apply(gDistance(pts_sp,z3,byid=TRUE),2,which.min)
cols = apply(gDistance(pts_sp, lakes, byid=TRUE), 2, which.min)
class(lakes)
class(pts_sp)
plot(pts_sp,pch=16,col=cols+1,add=TRUE)
plot(pts_sp,pch=1,add=TRUE)
with the work being done by gDistance(). The problem was that WKT wants
coordinates as "(123.456 789.012)" in parentheses separated by a space,
and the mailer line break (for my mailer, maybe nobody else saw this) hit
that critical space.
The main point was actually to say that gDistance() should also handle the
points to line question that came up today too:
gDistance(pts_sp, as(lakes, "SpatialLines"), byid=TRUE)
cols = apply(gDistance(pts_sp, as(lakes, "SpatialLines"), byid=TRUE), 2,
which.min)
plot(as(lakes, "SpatialLines"), lwd=2, col=1:5+1)
plot(pts_sp,pch=16,col=cols+1,add=TRUE)
plot(pts_sp,pch=1,add=TRUE)
Roger
>
--
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