On Wed, 22 Sep 2004, Johannes SCHNITZLER wrote:
> Dear All,
>
> In the moment i'm using the map and maptools package to read shapefiles
> and display the maps.
>
> I'm looking for the possibility to draw points (randomly positioned or
> positioned according to a grid) into the polygons instead of filling the
> polygons with colors.
>
I apologize for not replying earlier. I think that you can combine the
polylist class from the maptools package with the csr function in the
splancs package to position points randomly (the gridpts() finction can
be used instead of csr() to give grid points):
> library(maptools)
> library(splancs)
> source("dotsinpolys.R")
> x <- read.shape(system.file("shapes/sids.shp",
package="maptools")[1])
> ncpolys <- Map2poly(x, raw=FALSE)
> names(x$att.data)
> plot(ncpolys)
> try1 <- dotsInPolys(ncpolys, x$att.data$SID74)
> xx <- lapply(try1, function(x) {if (!is.null(x)) points(x, pch=18,
+ col="red")})> try2 <- dotsInPolys(ncpolys, x$att.data$SID74, f=gridpts)
> plot(ncpolys)
> xx <- lapply(try2, function(x) {if (!is.null(x)) points(x, pch=18,
+ col="red")})
where dotsinpolys.R is:
dotsInPolys <- function(pl, x, f=csr) {
if (!inherits(pl, "polylist")) stop("not a polylist
object")
if (length(pl) != length(x)) stop("different lengths")
n <- length(pl)
res <- vector(mode="list", length=n)
for (i in 1:n) {
if (x[i] > 0) {
if (attr(pl[[i]], "nParts") == 1) {
if (attr(pl[[i]], "ringDir") != 1)
warning(paste("hole with content at:", i))
res[[i]] <- f(matrix(c(ncpolys[[i]]), ncol=2), x[i])
} else {
areas <- rep(0, attr(pl[[i]], "nParts"))
for (j in 1:attr(pl[[i]], "nParts")) {
if (attr(pl[[i]], "ringDir")[j] == 1) {
from <- attr(pl[[i]], "pstart")$from[j]
to <- attr(pl[[i]], "pstart")$to[j]
areas[j] <- areapl(matrix(c(ncpolys[[i]][from:to,]),
ncol=2))
}
}
pareas <- areas/sum(areas)
px <- as.integer(round(pareas*x[i], digits = 0))
for (j in 1:attr(pl[[i]], "nParts")) {
if (px[j] > 0) {
from <- attr(pl[[i]], "pstart")$from[j]
to <- attr(pl[[i]], "pstart")$to[j]
pj <- matrix(c(ncpolys[[i]][from:to,]), ncol=2)
res[[i]] <- rbind(res[[i]], f(pj, px[j]))
}
}
}
res[[i]] <- matrix(res[[i]], ncol=2)
}
}
res
}
The extra work is because not all polygon units are a single polygon, and
some may be "lakes" within other polygons, doubling up on area. I
suspect
that the more sophisticated polygons in spatstat may be able to prevent
sampled points landing in the water. In multi-polygon units, points are
divided by relative area, and because of rounding may not be exactly the
same as the input count. The gridded points are also not likely to add up
exactly.
Could you indicate whether this function is any use, and if you would like
it added to some suitable package?
Roger Bivand
>
>
> For example:
>
> a map (shapefile) with 10 countries, 15 points in the polygon of country
> A, 20 points in the polygon of country B.....
>
>
>
>
>
> Thank you in advance for your help
>
>
>
> Johannes
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
>
--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Breiviksveien 40, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 93 93
e-mail: Roger.Bivand at nhh.no