Schnittbereiche in R extrahieren

19

Ich habe zwei Polygone. Eines enthält Felder (X, Y, Z) und das andere enthält Bodentypen (A, B, C, D). Ich möchte wissen, in welchem ​​Gebiet jedes Feldes sich welcher Bodentyp befindet. Ich habe folgendes versucht:

Bildbeschreibung hier eingeben

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

und gute Ergebnisse damit erzielt, mir mitzuteilen, welches Feld welchen Bodentyp enthält. Wie bekomme ich stattdessen die Fläche?

user2386786
quelle
1
Hinweis: st_intersection funktioniert nicht, wenn Ihre Punkte Längen- und Breitengrade sind. Sie haben keine geografischen Koordinaten angegeben, obwohl dies angedeutet ist, da es sich um Bodentypen handelt.
Fourier

Antworten:

24

Diese Methode verwendet die intersect()Funktion aus dem rasterPaket. Die Beispieldaten, die ich verwendet habe, sind nicht ideal (zum einen sind sie in nicht projizierten Koordinaten), aber ich denke, sie vermitteln die Idee.

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

Ergebnisse:

    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
quelle
2
Zur Verdeutlichung: Ich bevorzuge raster::intersectes, rgeos::gIntersectionweil ersteres die Attributdaten der beiden SpatialPolgonsDataFrameObjekte verknüpft, während letzteres die Attributdaten zu löschen scheint.
Matt SM
Vielen Dank für die vielen Details und die richtige Antwort. Sie haben mir sehr geholfen!!!
user2386786
4
Wenn Sie in "gIntersection" byid = TRUE verwenden, wird das Attribut IDS zurückgegeben, das mit merge zum Zuordnen der Attribute verwendet werden kann. Die Funktionen sind unterschiedlich und es sollte beachtet werden wie. Die Funktion "intersect" verwendet die überlappenden Ausdehnungen, während "gIntersection" die explizite Schnittmenge der Vektorgeometrien ist. Der Schnittpunkt-Ansatz ist ein quadratischer / rechteckiger Schnittpunkt und kein Schnittpunkt der tatsächlichen Polygone. Die Ausdehnung kann mit Ausdehnung und bbox neu definiert werden. Beide Ansätze bieten Vorteile.
Jeffrey Evans
1
@ JeffreyEvans Guter Punkt bezüglich gIntersection; Die Eingabe-Feature-IDs werden jedoch nicht direkt bereitgestellt, sondern verkettet und in der Feature-ID der Ausgabe gespeichert. Dies bedeutet die zusätzlichen Schritte zum Parsen der IDs und anschließenden Verknüpfen der Attribute. Ich möchte raster::intersectdiese Eingabe-IDs als zusätzliche Attribute in die Ausgabe aufnehmen.
Matt SM
1
Vielen Dank für den Hinweis, ich habe intersect_sp komplett verpasst. Interessanterweise verwendet es gIntersects. Schöne Abkürzung, wenn Sie möchten, dass die Attribute verbunden werden.
Jeffrey Evans
23

Hier ist ein alternativer Ansatz mit dem neuen sfPaket, das ersetzt werden soll sp. Alles ist viel sauberer und pfeifenfreundlicher:

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

Bildbeschreibung hier eingeben

   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
quelle