0

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)

Illustration of the overlapping polygons

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

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
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.