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