divide polygon shapefile into 3 equal areas

This probably on the limit of acceptable LOCs on this list but here goes:

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(paste("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)


if not convinced:

gArea(parts[[1]])

[1] 3375

gArea(parts[[2]])

[1] 3375.001

gArea(parts[[3]])

[1] 3374.999

Can easily chop into quarters too... There's some redundancy in the
code, and I'm sure it can be improved...

Barry




On Mon, Feb 29, 2016 at 6:14 PM, Boris Steipe <boris.steipe at utoronto.ca>
wrote:

Sounds like a fun little bit of code to write:

 - write a function that will return the area of a slice as a function

of a parameter X that can vary between some bounds on your shape: left to
right, or top to bottom etc. E.g. if you want to slice vertically, this
could be the area of the part of your polygon between the leftmost point
and a vertical line at X. (Adapt from here perhaps:
https://stat.ethz.ch/pipermail/r-sig-geo/2015-July/023168.html)

 - find the roots of that function for f(X, shape) - 1/3 * totalArea and

f(X, shape) - 2/3 * totalArea

   (

https://stat.ethz.ch/R-manual/R-devel/library/stats/html/uniroot.html )

B.

On Feb 29, 2016, at 12:57 PM, Shane Carey <careyshan at gmail.com> wrote:

ok thanks!!

I would like to slice it vertically and have 3 distinct areas of equal
area. So I need to chop it up into 3 areas of equal size essentially.

There is no tool to do it in QGIS!!

Thanks

On Mon, Feb 29, 2016 at 5:51 PM, Barry Rowlingson <
b.rowlingson at lancaster.ac.uk> wrote:

On Mon, Feb 29, 2016 at 5:37 PM, Shane Carey <careyshan at gmail.com>

wrote:

Hi,

Is it possible to divide a polygon into 3 equal areas using R?

Yes, in an infinite number of ways. Want to narrow it down?

Specifically, you could slice it vertically, horizontally, or at any
angle between. You could chop it into squares and reassign them (did
you want **contiguous** areas?). You could choose a point and three
radii angles that divide the polygon into 3 equal areas in an infinite
number of ways.

The rgeos package will help you chop polygons up, and then uniroot
can find the coordinates of lines or radii of angles that chop the
polygon first into 1/3 & 2/3 then chop the 2/3 into 1/2 and 1/2,
giving you three equal pieces.

I cant seem to be able to do it in QGIS.

If it can be done in R it can be done in Python and then it can be
done in QGIS...

Barry

--
Shane

      [[alternative HTML version deleted]]