# Splitting an administrative region into given number of parts

I have the following plot.

I want to divide the administrative region into given number (say 20) equal parts.

``````library(maps)
data(world.cities)

Pakistan <- data.frame(map("world", "Pakistan", plot = FALSE)[c("x","y")])

library(ggplot2)

p <- ggplot(Pakistan, aes(x = x, y = y)) +
geom_path(colour = 'green', linestyle = 2) +
coord_map() +
theme_bw() +
labs(x=" ", y=" ") +
theme(
panel.grid.minor = element_blank()
, panel.grid.major=element_blank()
, axis.ticks = element_blank()
, axis.text.x = element_blank()
, axis.text.y = element_blank()
, panel.border = element_blank()
)

print(p)
``````

Edited

I found this one

``````makeVchopper <- function(pol){
bb = bbox(pol)
delta = (bb[2,2] - bb[2,1])/10
xmin = bb[1,1]-delta
ymin = bb[2,1]-delta
ymax = bb[2,2]+delta

choppoly = function(xmax){
readWKT(sprintf("POLYGON((%s %s, %s %s, %s %s, %s %s, %s %s))",
xmin,ymin, xmin,ymax, xmax,ymax, xmax,ymin, xmin,ymin))
}
choppoly
}

slicer <- function(pol, xmin, xmax){
bb = bbox(pol)
delta = (bb[2,2] - bb[2,1])/10
ymax = bb[2,2] + delta
ymin = bb[2,1] - delta
r = readWKT(sprintf("POLYGON((%s %s, %s %s, %s %s, %s %s, %s %s))",
xmin,ymin, xmin,ymax, xmax,ymax, xmax,ymin, xmin,ymin))
gIntersection(pol,r)
}

chop_thirds <- function(pol, fractions=c(1/3, 2/3)){
chopper = makeVchopper(pol)
bb = bbox(pol)
xmin = bb[1,1]
xmax = bb[1,2]

totalArea = gArea(pol)

chopped_area = function(x){
gArea(gIntersection(chopper(x),pol))
}

edges = lapply(fractions, function(fraction){
target = totalArea * fraction
target_function = function(x){
chopped_area(x) - target
}
uniroot(target_function, lower=xmin, upper=xmax)\$root
})

xdelta = (xmax-xmin)/10
chops = matrix(c(xmin-xdelta, rep(edges,rep(2,length(edges))),
xmax+xdelta), ncol=2, byrow=TRUE)
apply(chops, 1, function(edges){
slicer(pol, edges[1], edges[2])
})

}

# Usage:

library(rgeos)
library(sp)

# sample data
pol <- readWKT("POLYGON((-180 -20, -140 55, 10 0, -140 -60, -180 -20), (-150 -20, -100 -10, -110 20, -150 -20))")
plot(pol)

# now split

parts = chop_thirds(pol)
plot(pol)

gArea(parts[[1]])
gArea(parts[[2]])
gArea(parts[[3]])
``````
• That looks like my code. Does it work? What's your problem? – Spacedman May 20 '19 at 14:45
• Thanks @Spacedman for your code. Would you like to guide me how to divide my shapefile into any number of parts. Thanks – MYaseen208 May 20 '19 at 15:01
• Your `Pakistan` object is made up of several line segments separated by `NA` values at rows 8, 243, 317, 384 and 542. You first need to make a continuous polygonal structure. – Spacedman May 20 '19 at 17:14
• Thanks @Spacedman for your comment. Wonder how to do what you are suggesting. – MYaseen208 May 20 '19 at 17:16
• Can you get a shapefile of Pakistan? Can you use the GADM data, like `p = raster::getData("GADM",country="PAK",level=0)`? – Spacedman May 20 '19 at 17:20

Working with an outline of Pakistan from the GADM data:

``````p = raster::getData("GADM",country="PAK",level=0)
``````

`p` is a SpatialPolygons data frame in lat-long coordinates. You should convert this to a proper equal-area projection but EPSG:3857 should be close enough:

``````pp = sp::spTransform(p, "+init=epsg:3857")
``````

Now chop using the code you pasted in your question. Specify a sequence of fractions from 1/20 to 19/20 get the split you want. This may take a while (twenty minutes?):

``````parts = chop_thirds(pp, fractions=seq(1/20, 19/20, len=19))
``````

Now plot the parts in 20 fabulous rainbow colours:

``````plot(pp)
cols = rainbow(20)
``````> range(unlist(lapply(parts, gArea)))