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
并取得了很好的结果，告诉我哪个字段包含哪种土壤类型。但是，我该如何获取该区域？
#1 楼
此方法使用intersect()软件包中的raster函数。我使用的示例数据并不理想（一方面，它们处于未投影的坐标中），但我认为它可以使您理解。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)
结果：
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
评论
需要说明的是：与rgeos :: gIntersection相比，我更喜欢raster :: intersect，因为前者将来自两个SpatialPolgonsDataFrame对象的属性数据联接在一起，而后者似乎会丢弃属性数据。
–马特SM
15年3月27日在1:03
感谢您提供许多详细信息和正确答案。你帮了我很多！！！
–user2386786
15年3月27日在8:01
如果在“ gIntersection”中使用byid = TRUE，它将返回属性IDS，该IDS可与merge一起使用以关联属性。功能不同，应注意如何操作。 “相交”功能使用重叠范围，而“ gIntersection”是矢量几何的显式交集。相交方法是正方形/矩形相交，而不是实际多边形的相交。可以使用扩展区和bbox重新定义扩展区。两种方法都有优点。
–杰弗里·埃文斯（Jeffrey Evans）
15年6月12日在18:22
@JeffreyEvans好一点的gIntersection;但是，输入要素ID不会直接提供，而是被串联并存储在输出的要素ID中。这意味着解析ID，然后加入属性的额外步骤。我确实希望raster :: intersect将这些输入ID作为附加属性包含在输出中。
–马特SM
15年6月12日在19:04
感谢您指出这一点，我完全错过了intersect_sp。有趣的是，它使用gIntersects。如果您想加入属性，则是不错的捷径。
–杰弗里·埃文斯（Jeffrey Evans）
2015年6月12日在20:46
#2 楼
这是使用新的sf软件包的替代方法，该软件包将替换sp。一切都更加清洁，而且对管道友好：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
评论
请注意，如果您的点是经度和纬度，则st_intersection将不起作用。您没有指定具有地理坐标，但是由于您在谈论土壤类型而被暗示。