rdataframespatialintersect

Intersect in R - miss one polygon

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!

Solution

  • 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.