Intersection d'une grille avec des cercles

Question

J'ai le problème suivant. J'ai des données d'un canton de CH (Fribourg) et j'aimerais croiser 3 cercles avec une cellule de grille créée dans cette zone. Voici un code pour générer la carte ci-jointe et le polygone:

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

Ce que je veux essentiellement à faire est d'obtenir la valeur 1 pour toutes les cellules de la grille en dehors des cercles, et 2 pour les cellules de la grille à l'intérieur des cercles. Ce qui le rend un peu compliqué, ce sont les cellules de la grille qui ne sont pas entièrement dans le cercle. Pour ces cellules, je voudrais attribuer la moyenne pondérée de 1 et 2 pondérée par la surface à l'extérieur et à l'intérieur du cercle respectivement.

Par exemple permet de dire que nous avons une cellule de grille en chevauchement avec un cercle comme l'image suivante (le rouge est la partie extérieure du cercle et bleu à l'intérieur)

Je veux que cette cellule pour prendre la valeur (1 * | A | + 2 * | B |)/2

Quelqu'un a-t-il une idée?

Answer 1

Vous pouvez utiliser la fonction rasterize dans la bibliothèque raster, avec le paramètre getCover=TRUE.

Si TRUE, la fraction de chaque cellule de la grille qui est couvert par les polygones est renvoyée (et les valeurs de champ, amusement, masque, et mise à jour sont ignorés. La fraction visée est estimée en divisant chaque cellule dans 100 sous-cellules et déterminer la présence/absence du polygone dans le centre de chaque sous-cellule

Si vous voulez quelque chose de plus précis, vous pouvez vous disaggregate la trame dans une trame plus précise (supérieure à 100, sinon il est le même que getCover), et récupérez la correspondance entre votre raster d'origine et le raster.