[R] divide polygon shapefile into 3 equal areas

[Shane Carey]

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

[Bert Gunter]

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 )

[Barry Rowlingson]

On Mon, Feb 29, 2016 at 5:37 PM, Shane Carey <care...@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 Carey]

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.rowl...@lancaster.ac.uk> wrote:

> On Mon, Feb 29, 2016 at 5:37 PM, Shane Carey <care...@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]]

[Boris Steipe]

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.

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

[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