27

I have two polygons. One contains fields(X,Y,Z) and the other contains soil types (A,B,C,D). I want to know what area of every field contains which type of soil. I tried the following:

enter image description here

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

and achieved good results with it telling me which field contains which soil type. However, how do I get the area instead?

user2386786
  • 549
  • 1
  • 5
  • 12
  • 2
    As a note, st_intersection won't work if your points are latitude and longitudes. You didn't specify that you had geographic coordinates, though it's hinted at since you are talking about soil types. – Fourier Mar 09 '19 at 19:27

2 Answers2

37

Here's an alternate approach using the new sf package, which is meant to replace sp. Everything is much cleaner, and pipe friendly:

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

enter image description here

   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
chenga
  • 3
  • 1
Matt SM
  • 1,841
  • 15
  • 23
29

This method uses the intersect() function from the raster package. The example data I've used aren't ideal (for one thing they're in unprojected coordinates), but I think it gets the idea across.

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

Results:

    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
Robert Hijmans
  • 10,683
  • 25
  • 35
Matt SM
  • 1,841
  • 15
  • 23
  • 2
    To clarify: I prefer raster::intersect over rgeos::gIntersection because the former joins the attribute data from the two SpatialPolgonsDataFrame objects, while the latter seems to drop the attribute data. – Matt SM Mar 27 '15 at 01:03
  • Thanks for the many details and the correct answer. You helped me a lot!!! – user2386786 Mar 27 '15 at 08:01
  • 5
    If you use byid=TRUE in "gIntersection" it will return attribute IDS which can be used with merge to associate the attributes. The functions are different and it should be noted how. The "intersect" function uses the overlapping extents whereas, "gIntersection" is the explicit intersection of the vector geometries. The intersect approach is a square/rectangular intersection and not an intersection of the actual polygons. The extent can be redefined using extent and bbox. There are advantages to both approaches. – Jeffrey Evans Jun 12 '15 at 18:22
  • 1
    @JeffreyEvans Good point re gIntersection; however, the input feature IDs aren't directly provided, they're concatenated and stored in the feature ID of the output. This means the extra steps of parsing the IDs, then joining in the attributes. I do wish raster::intersect included these input IDs as additional attributes in the output. – Matt SM Jun 12 '15 at 19:04
  • @JeffreyEvans Also, perhaps my interpretation of your comments re raster::intersect is incorrect; however, I disagree that it only uses overlapping extents. In fact, it does the true intersection of the vector geometries just as gIntersection. – Matt SM Jun 12 '15 at 19:06
  • @Matt SM, from the raster intersect help "Extent objects: Returns the intersection, i.e. the area of overlap of two Extent objects" and "SpatialPolygons* objects: Only the overlapping areas (if any) are returned". I also looked at the source and it is using the extent class object to define the spatial subset. – Jeffrey Evans Jun 12 '15 at 19:12
  • @JeffreyEvans Again, maybe I'm misunderstanding, but getMethod(intersect, c("SpatialPolygons", "SpatialPolygons")) seems contrary to what you're saying. Also, continuing from the example in my answer above, the two methods return same geometry, only difference is attributes: pi <- intersect(soil, field); plot(pi); pg <- gIntersection(soil, field, byid=T); plot(pg) – Matt SM Jun 12 '15 at 19:31
  • I see what you're saying re the source code, intersect.R uses extents only, but intersect_sp.R uses the full geometry. – Matt SM Jun 12 '15 at 19:34
  • 1
    Thanks for pointing that out, I completely missed intersect_sp. Interestingly, it uses gIntersects. Nice short cut if you want the attributes joined. – Jeffrey Evans Jun 12 '15 at 20:46
  • This code works fine without loading the maptools library (explicitly), at least on my system. – Daan Apr 30 '21 at 12:26