Re: [R] divide polygon shapefile into 3 equal areas
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.rowling...@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> 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 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 > 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 > >>
Re: [R] 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 Steipewrote: > 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 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 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
Re: [R] divide polygon shapefile into 3 equal areas
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 Careywrote: > 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 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.
Re: [R] divide polygon shapefile into 3 equal areas
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 Careywrote: > > 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.
Re: [R] divide polygon shapefile into 3 equal areas
On Mon, Feb 29, 2016 at 5:37 PM, Shane Careywrote: > 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 __ 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.
Re: [R] divide polygon shapefile into 3 equal areas
The r-sig-geo list might be a better place for this post. Cheers, Bert Bert Gunter "The trouble with having an open mind is that people keep coming along and sticking things into it." -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip ) On Mon, Feb 29, 2016 at 9:37 AM, Shane Careywrote: > Hi, > > Is it possible to divide a polygon into 3 equal areas using R? > > I cant seem to be able to do it in QGIS. > > Thanks > > -- > 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.