[R] [FORGED] intersection of two polygons which are not shapefiles

[Rolf Turner]

On 19/07/16 01:16, Adrienne Wootten wrote:
> All,
>
> Greetings! I hope things are going well for all!  I apologize if someone's
> already answered this and I'm just not finding it.  Here's what I'd like to
> do, but I'm hitting a brick wall with it.
>
> I have two sets of points for which I've already determined which ones
> points for the boundaries with the chull function.  What I need for what
> I'm doing is the coordinates where the two resulting polygons overlap.
> These are not raster grids, nor shapefiles.  What I have to work with are
> the two data frames with the points at the vertices of each polygon
> (included below).
>
>> chx1
>               x          y
> 1  0.5822569 -0.5555878
> 2  0.5338428 -0.5883604
> 3 -0.3442943 -0.6701115
> 4  -0.7409293  0.2286962
> 5  0.2147221  0.8061485
> 6  0.4914146  0.4941865
> 7  0.5822569 -0.5555878
>
>> chx2
>               x          y
> 1  0.7163506 -0.4357497
> 2  0.6513128 -0.5395180
> 3   0.1818315 -0.6317423
> 4  -0.6394281 -0.5765610
> 5 -0.6044681  0.1831627
> 6 -0.5799485  0.3231473
> 7  0.2248476  0.9601908
> 8  0.7163506 -0.4357497
>
>
> If anyone has any ideas on how to get what I'm after I'd appreciate it!
> I've tried a lot of things from the raster, rgeos, and more.  Knowing me
> it's something obvious I'm just not seeing right now.

You can do this very easily using the spatstat package:

library(spatstat)
X1 <- read.table(textConnection("
              x          y
1  0.5822569 -0.5555878
2  0.5338428 -0.5883604
3 -0.3442943 -0.6701115
4  -0.7409293  0.2286962
5  0.2147221  0.8061485
6  0.4914146  0.4941865
7  0.5822569 -0.5555878"))

X2 <- read.table(textConnection("
                x          y
1  0.7163506 -0.4357497
2  0.6513128 -0.5395180
3   0.1818315 -0.6317423
4  -0.6394281 -0.5765610
5 -0.6044681  0.1831627
6 -0.5799485  0.3231473
7  0.2248476  0.9601908
8  0.7163506 -0.4357497"))

X1 <- reverse.xypolygon(X1) # Vertices are in the wrong
X2 <- reverse.xypolygon(X2) # (clockwise) order.
W1 <- owin(poly=X1)
W2 <- owin(poly=X2)
WI <- intersect.owin(W1,W2)

plot(union.owin(W1,W2),col="blue",main="")
plot(WI,add=TRUE,col="red")

HTH

cheers,

Rolf Turner

P.S.  To extract the coordinates of the intersection polygon you can do:

WI$bdry[[1]]

R. T.

-- 
Technical Editor ANZJS
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276