在 R 中相交 - 错过一个多边形

Question

1.问题

我正在尝试在 R 中提取两个多边形形状的交集。第一个是分水岭多边形“ws_polygon_2”，第二个是 5 个雨量计的 voronoi 多边形，由 excel 表“DATA.xlsx”构建"，两者都可以在这里找到：link。

代码如下：

#[1] Montagem da tabela de coordenadas dos postos pluviométricos
library(sp)
library(readxl)
dados_precipitacao_1985 <- read_excel(path="C:/Users/.../DATA.xlsx") 
coordinates(dados_precipitacao_1985) <- ~ x + y 
proj4string(dados_precipitacao_1985) <- CRS("+proj=longlat +datum=WGS84") 
d_prec <- spTransform(dados_precipitacao_1985, CRSobj = "+init=epsg:3857") 

#[2] Coleta dos dados espaciais da bacia hidrográfica
library(rgdal)
bacia_Caio_Prado <- readOGR(dsn="C:/Users/...", layer="ws_polygon_2")
bacia_WGS <- spTransform(bacia_Caio_Prado, CRSobj = "+proj=longlat +datum=WGS84")
bacia_UTM <- spTransform(bacia_Caio_Prado, CRSobj = "+init=epsg:3857")

#[3] Poligonos de Thiessen - 1 INTERPOLAÇÃO
library(dismo)
library(rgeos)
library(raster)
library(mapview)
limits_voronoi_WGS <- c(-40.00,-38.90,-5.00,-4.50)
v_WGS <- voronoi(dados_precipitacao_1985, ext=limits_voronoi_WGS)
bc <- aggregate(bacia_WGS)
u_WGS_1 <- gIntersection(spgeom1 = v_WGS, spgeom2 = bc,byid=TRUE)
u_WGS_2 <- intersect(bc, v_WGS)

当我应用intersect 函数时，返回的变量u_WGS_2 是一个空间多边形数据框，只有4 个特征，而不是5 个。voronoi 对象v_WGS 也有5 个特征。

另一方面，当我应用 gIntesection 函数时，我得到了 5 个特性。但是，u_WGS_1 对象只是一个空间多边形，我丢失了降雨数据。

我想知道我是否犯了任何错误，或者是否有任何方法可以通过intersect 函数将 5 个要素与空间多边形数据框中的降雨数据聚合在一起。

我的目标是稍后通过rasterize 函数将这个空间多边形数据框与栅格中每个 voronoi 多边形的降雨数据进行转换，以便与其他插值结果和卫星数据进行比较。

看看这些结果。第一个是当我得到我想要的 SPDF（空间多边形数据框），但缺少 5º 功能。第二个是我获得的所有我想要的功能，但缺少降雨数据。 spplot(u_WGS_2, 'JAN') plot(u_WGS_1)

2。我尝试过的

1. 我查看ws_polygon_2 形状，寻找任何其他会污染形状并引导此结果的不需要的多边形。 形状仅由一个多边形特征组成，即正确的分水岭特征。

2. 我尝试使用aggregate 函数，如上所述，正如我在this 教程中看到的那样。 但我得到了相同的结果。

3. 我尝试使用 de u_WGS_1 和 d_precSpatial Point Data Frame 对象创建 SPDF。 实际上，我正在努力。如果这是解决我问题的正确答案，请帮我写一些代码。

谢谢！

Answer 1

当使用来自 sf 的st_intersection() 时，这不是问题，它保留了两个数据集中的数据。请注意，dismo::voronoi() 仅与 sp 对象兼容，因此降水数据需要以该格式提供，至少是暂时的。如果您对 sf 不满意，并且希望在实际交集之后继续使用 Spatial* 对象，只需在输出 sf 对象上调用 as() 方法，如图所示以下。

library(sf)

#[1] Montagem da tabela de coordenadas dos postos pluviométricos
dados_precipitacao_1985 <- readxl::read_excel(path="data/DATA.xlsx") 
dados_precipitacao_1985 <- st_as_sf(dados_precipitacao_1985, coords = c("x", "y"), crs = 4326)
dados_precipitacao_1985_sp <- as(dados_precipitacao_1985, "Spatial")

#[2] Coleta dos dados espaciais da bacia hidrográfica
bacia_Caio_Prado <- st_read(dsn="data/SHAPE_CORRIGIDO", layer="ws_polygon_2")

#[3] Poligonos de Thiessen - 1 INTERPOLAÇÃO
limits_voronoi_WGS <- c(-40.00,-38.90,-5.00,-4.50)
v_WGS <- dismo::voronoi(dados_precipitacao_1985_sp, ext=limits_voronoi_WGS)
v_WGS_sf <- st_as_sf(v_WGS)

u_WGS_3 <- st_intersection(bacia_Caio_Prado, v_WGS_sf)
plot(u_WGS_3[, 6], key.pos = 1)

Answer 2

缺失的多边形被移除，因为它无效

library(raster)
bacia <- shapefile("SHAPE_CORRIGIDO/ws_polygon_2.shp")
rgeos::gIsValid(bacia)

#[1] FALSE
#Warning message:
#In RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid") :
#  Ring Self-intersection at or near point -39.070555560000003 -4.8419444399999998

自交点在这里：

zoom(bacia, ext=extent(-39.07828, -39.06074, -4.85128, -4.83396))
points(cbind( -39.070555560000003, -4.8419444399999998))

无效的多边形被删除，因为它们被认为是由 intersect 生成的。在这种情况下，无效数据已经存在并且应该被保留。我会看看我能不能解决这个问题。