Intersect grid with circles

Question

I have the following problem. I have data from a canton of CH (Fribourg) and I would like to intersect 3 circles with a grid cell created in this area. Here is some code to generate the map attached and the polygon:

data file

library(raster)

centroid.x <- c(566052, 576768, 560652)
centroid.y <- c(155501, 185501, 165325)`

centroids <- data.frame(centroid.x = centroid.x, centroid.y = centroid.y)

radius <- 1000
n = 1000
pts = seq(0, 2 * pi, length.out = n)
ZH = cbind(centroid.x[1] + radius * sin(pts), centroid.y[1] + radius * cos(pts))
WH <- cbind(centroid.x[2] + radius * sin(pts), centroid.y[2] + radius * cos(pts))
GL <- cbind(centroid.x[3] + radius * sin(pts), centroid.y[3] + radius * cos(pts))
sl = SpatialPolygons(list(Polygons(list(Polygon(ZH)), 1), Polygons(list(Polygon(WH)), 2) 
                          Polygons(list(Polygon(GL)), 3))) `

CRS_CH <- CRS("+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs")
crs(sl) <- CRS_CH
GRID <- 1000
grid <- raster(extent(pol))

res(grid) <- GRID
proj4string(grid)<-proj4string(pol)

gridpolygon <- rasterToPolygons(grid)
plot(gridpolygon)

dry.grid1000 <- intersect(pol, gridpolygon)
plot(dry.grid1000)
plot(sl, add = TRUE, border = "red")

What I basically want to do is get the value 1 for all grid cells outside the circles, and 2 for the grid cells inside the circles. What makes it a bit tricky is the grid cells that are not entirely in the circle. For these grid cells I would like to assign the weighted mean of 1 and 2 weighted by the area outside and inside the circle respectively.

For example lets say that we have a grid cell overlapping with a circle as the following image (red is the part outside the circle and blue inside)

I want this cell to take the value (1*|A| + 2*|B|)/2

Anyone any idea?

Answer 1

You can use function rasterize in library raster, with parameter getCover=TRUE.

If TRUE, the fraction of each grid cell that is covered by the polygons is returned (and the values of field, fun, mask, and update are ignored. The fraction covered is estimated by dividing each cell into 100 subcells and determining presence/absence of the polygon in the center of each subcell

If you want something more precise, you can disaggregate yourself the raster into a more precise raster (higher than 100, otherwise it is the same as getCover), and retrieve the correspondence between your original raster and the disaggregated one.