R SF – Count Number of Polygons in Multipolygon Using R SF

polygonrsf

How can I count the number of polygons that are present in a multipolygon with the R package sf? Example: In the following image that represents a multipolygon I count 3 polygons. I want to be able to store this information as an attribute. Below I added the code that creates the multipolygon.

# Create Multipolygon: https://cran.r-project.org/web/packages/sf/vignettes/sf1.html
library(sf)
p1 <- rbind(c(0,0), c(1,0), c(3,2), c(2,4), c(1,4), c(0,0))
p2 <- rbind(c(1,1), c(1,2), c(2,2), c(1,1))
pol <-st_polygon(list(p1,p2))
p3 <- rbind(c(3,0), c(4,0), c(4,1), c(3,1), c(3,0))
p4 <- rbind(c(3.3,0.3), c(3.8,0.3), c(3.8,0.8), c(3.3,0.8), c(3.3,0.3))[5:1,]
p5 <- rbind(c(3,3), c(4,2), c(4,3), c(3,3))
mpol <- st_multipolygon(list(list(p1,p2), list(p3,p4), list(p5)))

mpol_geom = st_sfc(mpol, crs = 4326)   

mpol_attrib = data.frame(                                   
name = "MULTIPOLYGON"
)

mpol_sf = st_sf(mpol_attrib, geometry = mpol_geom)   

plot(mpol_sf)

Best Answer

EDIT -- correcting the answer using your sf object, I note that

length(st_geometry(mpol_sf)[[1]]) or length(mpol_sf$geometry[[1]])

should work.

Related Solutions

R – Count Overlapping Polygons in Single Shapefile

In R, you can used the sp package and over function to do this. I adapted an example data set and the solution from this post by Roger Bivand.

Set up the example data:

library(sp)
library(rgeos)
library(raster)
library(rworldmap)

box <- readWKT("POLYGON((-180 90, 180 90, 180 -90, -180 -90, -180 90))")
proj4string(box) <- CRS("+proj=cea +datum=WGS84")
set.seed(1)
pts <- spsample(box, n=2000, type="random")
pols <- gBuffer(pts, byid=TRUE, width=50) # create circle polys around each point
merge = sample(1:40, 100, replace = T) # create vector of 100 rand #s between 0-40 to merge pols on

Sp.df1 <- gUnionCascaded(pols, id = merge) # combine polygons with the same 'merge' value
# create SPDF using polygons and randomly assigning 1 or 2 to each in the @data df
Sp.df <- SpatialPolygonsDataFrame(Sp.df1,
                    data.frame(z = factor(sample(1:2, length(Sp.df1), replace = TRUE)),
                                                    row.names= unique(merge)))
Sp.df <- crop(Sp.df, box)
colors <- c(rgb(r=0, g=0, blue=220, alpha=50, max=255), rgb(r=220, g=0, b=0, alpha=50, max=255))

land <- getMap()

overlay.map <- spplot(Sp.df, zcol = "z", col.regions = colors,
                   col = NA, alpha = 0.5, breaks=c(0,1),
                   sp.layout = list("sp.polygons", land, fill = "transparent",
                                    col = "grey50"))
overlay.map

Then to get the count (or other statistics)...

# find the count of polygons below each grid cell
GT <- GridTopology(c(-179.5, -89.5), c(1, 1), c(360, 180))
SG <- SpatialGrid(GT)
proj4string(SG) <- CRS("+proj=cea +datum=WGS84")


o <- over(SG, Sp.df1, returnList=TRUE)
ct <- sapply(o, length)
summary(ct)
SGDF <- SpatialGridDataFrame(SG, data=data.frame(ct=ct))
spplot(SGDF, "ct", col.regions=bpy.colors(20))

[GIS] Count Polygons in Polygons

To find which parcel is located in which zone use the Intersection Tool under Vector-> Geoprocessing Tools with both layers as input layers.

The result is a layer, which contains every parcel (or part of it) which is inside the zone layer and the attribute table contains all attributes from both input layers.

Now you can either directly work with the intersection layer or split it with Vector-> Data Management Tools-> Split vector layer using the field which specifies the outer zone. Each of the new file now has directly the number of features which are located inside the corresponding zone.

Best Answer

Related Solutions

R – Count Overlapping Polygons in Single Shapefile

[GIS] Count Polygons in Polygons

Related Question