Пересечь в R - пропустить один многоугольник

Question

1.Проблема

Я пытаюсь извлечь пересечение двух форм многоугольников в R. Первый - это полигон водораздела "ws_polygon_2", а второй - это полигоны вороной из 5 дождемеров, которые были построеныиз листа Excel "DATA.xlsx", оба доступны здесь: ссылка .

Код следующий:

#[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 является только пространственным многоугольником, и я теряю данные об осадках.

Я хотел бы знать, допускаю ли я какую-либо ошибку или есть ли способ объединить 5 объектов сданные об осадках в фрейме пространственных полигонов с помощью функции intersect.

Моя цель - преобразовать этот фрейм данных пространственного многоугольника с данными об осадках для каждого многоугольника вороного в растре с помощью функции rasterize позже, чтобы сравнить с другими интерполяционными результатами и спутниковыми данными.

Посмотри эти результаты.Первый - когда я получаю SPDF (пространственный полигональный фрейм данных), который мне нужен, но отсутствует функция 5º.Второй - тот, который я получаю со всеми необходимыми функциями, но пропускаю данные об осадках.spplot(u_WGS_2, 'JAN') plot(u_WGS_1)

2.То, что я пробовал

1. Я изучаю фигуру ws_polygon_2 в поисках любого другого нежелательного многоугольника, который бы искажал форму и направлял к этим результатам. Форма состоит только из одного многоугольника, правильного элемента водораздела.

2. Я пытался использовать функцию aggregate, как указано выше, и как я увиделв этом учебнике. Но я получил тот же результат.

3. Я пытался создать SPDF с объектом u_WGS_1 и d_prec Spatial Point Data Frame. На самом деле, я работаю над этим.И если это правильный ответ на мою проблему, пожалуйста, помогите мне с кодом.

Спасибо!

Answer 1

Отсутствующий многоугольник удален, потому что он недействителен

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

Неверные многоугольники удалены, так как предполагается, что они были полученыпересекаются.В этом случае недействительные данные уже были там и должны были быть сохранены.Я посмотрю, смогу ли я это исправить.

Answer 2

Это не проблема при использовании st_intersection() из sf , который сохраняет данные из обоих наборов данных.Помните, что dismo::voronoi() совместим только с объектами sp , поэтому данные об осадках должны быть доступны в этом формате, по крайней мере, временно.Если вы не чувствуете себя комфортно с sf и предпочитаете продолжать работу с объектами Spatial * после фактического пересечения, просто вызовите метод as() для выходного объекта sf , как показано ниже.

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)