[R] divide polygon shapefile into 3 equal areas

[Shane Carey]

This is super, thanks!! However, I cannot read in my shapefile. I am using
readShapeSpatial and the error I receive is: Error in getinfo.shape(fn):
Error opening SHP file. The projection is Irish Transverse Mercator!!


Thanks in advance

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

> 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.
>



-- 
Shane

	[[alternative HTML version deleted]]