Splitting overlapping polygons in a layer and creating new ones with merged attributes using R

Question

I have a large shapefile with plenty of overlapping polygons. From the areas where they intersect, I would like to create new polygons (and thus split the original ones). I would like to also to merge the stored attributes - please note here I have field_ID numeric ID that would need to be stored either as a character split by comma or as a vector. This should result in 6 polygons.

I am using R.

library(sp)
# Funciton creating a spatial polygon:
create_poly <- function(coordin,id_number){     
  coordinates = matrix(c(coordin[1], coordin[3],
                         coordin[1], coordin[4],
                         coordin[2], coordin[4],
                         coordin[2], coordin[3],
                         coordin[1], coordin[3]), 
                       ncol = 2, byrow = TRUE)
  P_1 = Polygon(coordinates)
  Ps_1 = SpatialPolygons(list(Polygons(list(P_1), ID = id_number)), proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
  ID_1 <- sapply(slot(Ps_1, "polygons"), function(x) slot(x, "ID"))
  df <- data.frame(rep(0, length(ID_1)), row.names=ID_1)
  Ps1_df <- SpatialPolygonsDataFrame(Ps_1, df)
  Ps1_df
}

# Create shapefile with polygon boundaries of interest:
Pol_1 <- create_poly(c(80,100,20,40),1)
Pol_2 <- create_poly(c(90,120,30,50),1)
Pol_3 <- create_poly(c(130,150,30,50),1)
Pol_4 <- create_poly(c(145,170,40,60),1)

Pol_1@data$properties <- c("red")
Pol_2@data$properties <- c("blue")
Pol_3@data$properties <- c("green")
Pol_4@data$properties <- c("purple")

Pol_1@data$field_ID <- 1
Pol_2@data$field_ID <- 2
Pol_3@data$field_ID <- 3
Pol_4@data$field_ID <- 4

# Join the polygons into one layer:
Pol_large <- do.call(rbind, list(Pol_1, Pol_2, Pol_3, Pol_4))

plot(Pol_large, axes = TRUE)

I can intersect those polygons but in result I do get 8 - the intersected areas are double (as they have two different properties) and the original polygons are not cut.

library(raster)
inter1=raster::intersect(Pol_large, Pol_large)

I also tried changing the polygons into lines and polygonising again, but it does loose the attributes columns and would need additional processing:

library(rgeos)
pol_merge = gLineMerge(as(Pol_large,"SpatialLines"))
regions = gPolygonize(gNode(pol_merge))
#Show the 6 regions:
par(mar=c(0,0,0,0))
par(mfrow=c(2,3))
for(i in 1:length(regions)){plot(regions);plot(regions[i,],add=TRUE,col="red")}

gContains(Pol_large, regions, byid=TRUE)

Answer 1

You can convert your Pol_large into an sf object using the library sf. And then apply the function st_intersection to it.

Pol_inters <- st_intersection(st_as_sf(Pol_large))

Your attribute table then looks like this, with the IDs separated by comma:

Simple feature collection with 6 features and 5 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 80 ymin: 20 xmax: 170 ymax: 60
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
  rep.0..length.ID_1.. properties field_ID                       geometry n.overlaps origins
1                    0        red        1 POLYGON ((80 20, 80 40, 90 ...          1       1
2                    0        red        1 POLYGON ((90 40, 100 40, 10...          2    1, 2
3                    0       blue        2 POLYGON ((90 40, 90 50, 120...          1       2
4                    0      green        3 POLYGON ((130 30, 130 50, 1...          1       3
5                    0      green        3 POLYGON ((145 50, 150 50, 1...          2    3, 4
6                    0     purple        4 POLYGON ((145 50, 145 60, 1...          1       4