CODEWITHSUNDEEP
rgridgisrasterr-raster
Akis
Akis

Reputation: 130

Intersect grid with circles

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")

enter image description here

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)

enter image description here

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

Anyone any idea?

Upvotes: 1

Views: 465

Answers (1)

S&#233;bastien Rochette
S&#233;bastien Rochette

Reputation: 6661

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.

Upvotes: 1

Related Questions