Intersect in R - miss one polygon

Question

1. The problem

I'm trying to extract the intersection of two polygons shapes in R. The first is the watershed polygon "ws_polygon_2", and the second is the Voronoi polygons of 5 rain gauges which was constructed from the Excel sheet "DATA.xlsx", both available here: link.

The code is the following:

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

When I apply the intersect function, the variable returned u_WGS_2 is a spatial polygon data frame with only 4 features, instead of 5. The Voronoi object v_WGS has 5 features as well.

By other hand, when I apply the gIntesection function, I get 5 features. However, the u_WGS_1 object is a spatial polygon only and I loss the rainfall data.

I'd like to know if I am committing any mistake or if there is any way to get the 5 features aggregated with the rainfall data in a spatial polygon data frame through the intersect function.

My objective is to transform this spatial polygon data frame with the rainfall data for each Voronoi polygon in a raster through the rasterize function later to compare with other interpolating results and satellite data.

Look these results. The first one is when I get the SPDF (Spatial Polygon Data Frame) I want, but missing the 5º feature. The second is the one I get with all the features I want, but missing the rainfall data. spplot(u_WGS_2, 'JAN') plot(u_WGS_1)

2. What I've tried

1. I look into the ws_polygon_2 shape searching for any other unwanted polygon who would pollute the shape and guide to this results. The shape is composed by only one polygon feature, the correct watershed feature.

2. I tried to use the aggregate function, as above, and as I saw in this tutorial. But I got the same result.

3. I tried to create a SPDF with de u_WGS_1 and the d_precSpatial Point Data Frame object. Actually, I'm working on it. And if it is the correct answer to my trouble, please help me with some code.

Thank you!

Answer 1

This is not an issue when using st_intersection() from sf, which retains the data from both data sets. Mind that dismo::voronoi() is compatible with sp objects only, so the precipitation data needs to be available in that format, at least temporarily. If you do not feel comfortable with sf and prefer to continue working with Spatial* objects after the actual intersection, simply invoke the as() method upon the output sf object as shown below.

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

The missing polygon is removed because it is invalid

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

The self-intersection is here:

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

Invalid polygons are removed as they are assumed to have been produced by intersect. In this case, the invalid data was already there and should have been retained. I will see if I can fix that.