Crear un círculo de radio determinado alrededor de un punto y, a continuación, encontrar el área de superposición en un shapefile

Question

Mi caso de uso está relacionado con el cálculo del área de un polígono en un archivo de forma que es cubierto por otro polígono.

Específicamente, quiero hacer lo siguiente:

1. Crear círculos que tienen un radio de 100 metros de radio, desde una base de datos.marco de long/lat puntos
2. La superposición de estos círculos en un shapefile
3. Calcular el área de cada uno de los polígonos que están cubiertos por estos círculos.
4. Idealmente volver un dataframe donde cada fila contiene el área cubierta por todos los círculos para cada polígono.

Aquí creo que algunos datos simulados puntos.

mean_lat <- 46.07998
sd_lat <- 0.1609196
mean_long <- 8.931849
sd_long <-  0.1024659

dat_sim <- data.frame(lat = rnorm(500, mean = mean_lat, sd = sd_lat),
                      long = rnorm(500, mean = mean_long, sd = sd_long))

A continuación se obtienen los shapefiles

library(raster)
library(tidyverse)
library(sf)

# regular sp format
dat_ticino_sp <- getData(name = "GADM", 
                         country = "CHE", 
                         level = 3) %>% 
  subset(.@data$NAME_1 == "Ticino")


# convert to sf format
dat_ticino_sf <- getData(name = "GADM", 
                         country = "CHE", 
                         level = 3) %>%
  # convert to simple features
  sf::st_as_sf() %>%
  # Filter down to Ticino
  dplyr::filter(NAME_1 == "Ticino")

Así que lo que quiero hacer es :

• crear círculos con un radio de 100 metros de radio en la dat_sim, para superponer sobre el shapefile dat_ticino_sf o dat_ticino_sp archivos.

• Para cada polígono en el shapefile, calcular cuánto de ese polígono es cubierto por los círculos.

• devolver un dataframe donde cada fila (polígono) contiene el área cubierta por los círculos.

• la parcela esta

Los puntos de bonificación si esto se puede hacer de Características Simples

Answer 1

Creo que debería hacer lo que quiera.

Básicamente lo que he dicho en mi comentario. Una cosa que no he mencionado es que se desea transformar a una igualdad de área de proyección que utiliza metros, así que usted puede búfer de 100m y estar seguro de que las áreas calculadas a través de su estudio son equivalentes. No sé la correcta para Suiza, pero me he encontrado con este MODO de respuesta y usado.

La belleza de la sf es que integra (casi por arte de magia!) tan bien con dplyr - ver el último paso.

library(raster)
library(tidyverse)
library(sf)
#> Linking to GEOS 3.5.0, GDAL 2.1.1, proj.4 4.9.3

# convert to sf format
dat_ticino_sf <- getData(name = "GADM", 
                         country = "CHE", 
                         level = 3) %>%
  # convert to simple features
  sf::st_as_sf() %>%
  # Filter down to Ticino
  dplyr::filter(NAME_1 == "Ticino") %>% 
  st_transform(3035) # Reproject to EPSG:3035 as mentioned above

mean_lat <- 46.07998
sd_lat <- 0.1609196
mean_long <- 8.931849
sd_long <-  0.1024659

set.seed(42)
dat_sim <- data.frame(lat = rnorm(500, mean = mean_lat, sd = sd_lat),
                      long = rnorm(500, mean = mean_long, sd = sd_long))

# Convert to sf, set the crs to EPSG:4326 (lat/long), 
# and transform to EPSG:3035
dat_sf <- st_as_sf(dat_sim, coords = c("long", "lat"), crs = 4326) %>% 
  st_transform(3035)

# Buffer circles by 100m
dat_circles <- st_buffer(dat_sf, dist = 100)

# Intersect the circles with the polygons
ticino_int_circles <- st_intersection(dat_ticino_sf, dat_circles)
#> Warning: attribute variables are assumed to be spatially constant
#> throughout all geometries

## Plot a zoomed in map of polygons and overlay intersected circles to double check:
## (I assumed NAME_3 is the best unique identifier for the polygons)
bb <- st_bbox(ticino_int_circles)

plot(dat_ticino_sf[, "NAME_3"], 
     xlim = c(mean(c(bb["xmin"], bb["xmax"])) - 1000, 
              mean(c(bb["xmin"], bb["xmax"])) + 1000), 
     ylim = c(mean(c(bb["ymin"], bb["ymax"])) - 1000, 
              mean(c(bb["ymin"], bb["ymax"])) + 1000))

plot(ticino_int_circles[, "NAME_3"], add = TRUE)

# Summarize by NAME_3
# which aggregates all of the circles by NAME_3, 
# then calculate the area
ticino_int_circles_summary <- group_by(ticino_int_circles, NAME_3) %>% 
  summarise() %>% 
  mutate(area = st_area(.))