1

I have the following plot.

enter image description here

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 enter image description here

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)
plot(parts[[1]], add=TRUE, col=1)
plot(parts[[2]], add=TRUE, col=2)
plot(parts[[3]], add=TRUE, col=3) 

gArea(parts[[1]])
gArea(parts[[2]])
gArea(parts[[3]])
  • 1
    That looks like my code. Does it work? What's your problem? – Spacedman May 20 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 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 at 17:14
  • Thanks @Spacedman for your comment. Wonder how to do what you are suggesting. – MYaseen208 May 20 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 at 17:20
2

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)
for(i in 1:20){plot(parts[[i]],col=cols[i],add=TRUE)}

enter image description here

Note that some of the parts are divided into two sections because the boundary is concave. You can see this at the farthest east and west and also three of the blue ones in the middle. But all the sums of the areas in each of the 20 colours are equal to within quite a high tolerance (which can be tuned in the code):

> range(unlist(lapply(parts, gArea)))
[1] 58698950917 58698950962
  • Thanks @Spacedman for very useful answer. I wonder if there is a method to divide administrative region into parts which take the respect of roads. – MYaseen208 May 21 at 3:18
  • If you can define precisely what you mean mathematically, topologically, and geometrically, then you could ask that as a question. – Spacedman May 21 at 11:22

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.