Extracción de áreas de intersección en R

Tengo dos polígonos Uno contiene campos (X, Y, Z) y el otro contiene tipos de suelo (A, B, C, D). Quiero saber qué área de cada campo contiene qué tipo de suelo. Intenté lo siguiente:

ingrese la descripción de la imagen aquí

library(rgdal)
library(rgeos)
Field<-readOGR("./","Field")
Soil<-readOGR("./","Soil")
Results<-gIntersects(Soil,Field,byid=TRUE)
rownames(Results)<-Field@data$FieldName
colnames(Results)<-Soil@data$SoilType

> Results
      A     B     C     D
Z  TRUE FALSE FALSE FALSE
Y FALSE  TRUE  TRUE FALSE
X  TRUE  TRUE  TRUE  TRUE

y logré buenos resultados diciéndome qué campo contiene qué tipo de suelo. Sin embargo, ¿cómo obtengo el área en su lugar?

shapefile polygon r intersection rgdal usuario2386786
fuente

Como nota, st_intersection no funcionará si sus puntos son latitud y longitud. No especificó que tenía coordenadas geográficas, aunque se insinúa ya que está hablando de tipos de suelo.

Fourier

Respuestas:

Este método usa la intersect()función del rasterpaquete. Los datos de ejemplo que he usado no son ideales (por un lado, están en coordenadas no proyectadas), pero creo que hace que la idea se transmita.

library(sp)
library(raster)
library(rgdal)
library(rgeos)
library(maptools)

# Example data from raster package
p1 <- shapefile(system.file("external/lux.shp", package="raster"))
# Remove attribute data
p1 <- as(p1, 'SpatialPolygons')
# Add in some fake soil type data
soil <- SpatialPolygonsDataFrame(p1, data.frame(soil=LETTERS[1:12]), match.ID=F)

# Field polygons
p2 <- union(as(extent(6, 6.4, 49.75, 50), 'SpatialPolygons'),
             as(extent(5.8, 6.2, 49.5, 49.7), 'SpatialPolygons'))
field <- SpatialPolygonsDataFrame(p2, data.frame(field=c('x','y')), match.ID=F)
projection(field) <- projection(soil)

# intersect from raster package
pi <- intersect(soil, field)
plot(soil, axes=T); plot(field, add=T); plot(pi, add=T, col='red')

# Extract areas from polygon objects then attach as attribute
pi$area <- area(pi) / 1000000

# For each field, get area per soil type
aggregate(area~field + soil, data=pi, FUN=sum)

Imgur

Resultados:

    field soil         area
1      x    A 2.457226e+01
2      x    B 2.095659e+02
3      x    C 5.714943e+00
4      y    C 5.311882e-03
5      x    D 7.620041e+01
6      x    E 3.101547e+01
7      x    F 1.019455e+02
8      x    H 7.106824e-03
9      y    H 2.973232e+00
10     y    I 1.752702e+02
11     y    J 1.886562e+02
12     y    K 1.538229e+02
13     x    L 1.321748e+02
14     y    L 1.182670e+01

Matt SM
fuente

Para aclarar: prefiero raster::intersecthacerlo rgeos::gIntersectionporque el primero se une a los datos de los atributos de los dos SpatialPolgonsDataFrameobjetos, mientras que el segundo parece descartar los datos de los atributos.

Matt SM

Gracias por los muchos detalles y la respuesta correcta. ¡¡¡Me ayudaste mucho!!!

user2386786

Si usa byid = TRUE en "gIntersection", devolverá IDS de atributo que se puede usar con merge para asociar los atributos. Las funciones son diferentes y debe tenerse en cuenta cómo. La función "intersectar" utiliza las extensiones superpuestas, mientras que "gIntersection" es la intersección explícita de las geometrías vectoriales. El enfoque de intersección es una intersección cuadrada / rectangular y no una intersección de los polígonos reales. La extensión se puede redefinir usando la extensión y bbox. Hay ventajas para ambos enfoques.

Jeffrey Evans

@JeffreyEvans Buen punto re gIntersection; sin embargo, las ID de funciones de entrada no se proporcionan directamente, se concatenan y almacenan en la ID de función de la salida. Esto significa los pasos adicionales de analizar las ID y luego unir los atributos. Desearía raster::intersectincluir estas ID de entrada como atributos adicionales en la salida.

Matt SM

Gracias por señalar eso, me perdí por completo intersect_sp. Curiosamente, utiliza gIntersects. Buen atajo si quieres unir los atributos.

Jeffrey Evans

Aquí hay un enfoque alternativo que usa el nuevo sfpaquete, que está destinado a reemplazar sp. Todo es mucho más limpio y fácil de usar:

library(sf)
library(tidyverse)

# example data from raster package
soil <- st_read(system.file("external/lux.shp", package="raster")) %>% 
  # add in some fake soil type data
  mutate(soil = LETTERS[c(1:6,1:6)]) %>% 
  select(soil)

# field polygons
field <- c("POLYGON((6 49.75,6 50,6.4 50,6.4 49.75,6 49.75))",
        "POLYGON((5.8 49.5,5.8 49.7,6.2 49.7,6.2 49.5,5.8 49.5))") %>% 
  st_as_sfc(crs = st_crs(soil)) %>% 
  st_sf(field = c('x','y'), geoms = ., stringsAsFactors = FALSE)

# intersect - note that sf is intelligent with attribute data!
pi <- st_intersection(soil, field)
plot(soil$geometry, axes = TRUE)
plot(field$geoms, add = TRUE)
plot(pi$geometry, add = TRUE, col = 'red')

# add in areas in m2
attArea <- pi %>% 
  mutate(area = st_area(.) %>% as.numeric())

# for each field, get area per soil type
attArea %>% 
  as_tibble() %>% 
  group_by(field, soil) %>% 
  summarize(area = sum(area))

   field  soil      area
   <chr> <chr>     <dbl>
1      x     A  24572264
2      x     B 209573036
3      x     C   5714943
4      x     D  76200409
5      x     E  31015469
6      x     F 234120314
7      y     B   2973232
8      y     C 175275520
9      y     D 188656204
10     y     E 153822938
11     y     F  11826698

Matt SM
fuente