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.ste...@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 <careys...@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.rowling...@lancaster.ac.uk> wrote:
>>
>>> On Mon, Feb 29, 2016 at 5:37 PM, Shane Carey <careys...@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@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@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@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.
