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:
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?
Answer
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)
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
No comments:
Post a Comment