How to overlay raster and point layers and calculate summary statistics

Question

I'm looking for some advice on basic geospatial statistics. I'm using a raster file from Worldpop , indicating the population in Brazil every 100m2. I have another latlong dataset with the coordinates of hospitals in Brazil.

I would like to do the following:

1. Identify areas (and possibly create polygons) of areas that are within 1km from the hospitals, between 2-10km and above 10km
2. Calculate the number of people in each of the zones above

I'd like to provide some reproducible example, but the raster file is very large. There are some instructions on how to do this with two separate lists of lat/lon points, but I can't figure out how to do this with a raster file.

Any ideas?

Answer 1

Example data

library(raster)
bra <- getData('GADM', country="BRA", level=1)
r <- raster(bra, res=1)
values(r) <- 1:ncell(r)
r <- mask(r, bra)
pts <- coordinates(bra)
# plot(r)
# points(pts)

Solution

b1 <- extract(r, pts, buffer=100000)  # 100 km
b2 <- extract(r, pts, buffer=200000)  # 200 km
pop1 <- sapply(b1, sum)
pop2 <- sapply(b2, function(i)sum(i, na.rm=TRUE)) - pop1

To see the areas

spts <- SpatialPoints(pts, proj4string=crs(bra))
buf1 <- buffer(spts, width=100000, dissolve=FALSE)
buf2 <- buffer(spts, width=200000, dissolve=FALSE)

# adding IDs so that they can also be used in "extract"
buf1 <- SpatialPolygonsDataFrame(buf1, data.frame(id1=1:length(buf1)))
buf2 <- SpatialPolygonsDataFrame(buf2, data.frame(id2=1:length(buf2)))

# To combine buf1 and buf2 you could do
# buf <- (buf2-buf1) + buf1
# but in this example there are overlapping buffers, so I do
bb <- list()
for (i in 1:length(buf1)) {
    bb[[i]] <- (buf2[i,]-buf1[i,]) + buf1[i,]
}
buf <- do.call(bind, bb)

plot(r)
plot(buf,  col=c("red", "blue"), add=TRUE)

And now you could do

z <- extract(r, buf, fun=sum, na.rm=TRUE)
z <- cbind(data.frame(buf), z)
head(z)

To get the same result as above for pop1 and pop2

head(pop1)
head(pop2)