Wednesday, 29 August 2018

shapefile - Extracting intersection areas in R


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?



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)

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

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...