[R] dot density maps

[Roger Bivand]

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