[R] divide polygon shapefile into 3 equal areas

[Barry Rowlingson]

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]]
>>
>> ______________________________________________
>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.