[R-sig-Geo] Create a mask for multiple polygons in spatstat package
Rolf Turner
r.turner at auckland.ac.nz
Fri Nov 1 23:05:35 CET 2013
Your email was admirably structured in terms of providing reproducible code
which was set out very clearly and tidily. However I am afraid that I
don't understand
what it is that you are trying to do.
I think perhaps that you are using incorrect or misleading terminology.
What do you
mean by saying that you are trying to create a *mask*? Do you mean that
you are
trying to create a *window*?
Note that union.owin() forms a *window* which is the union of two or
more window
objects (objects of class "owin").
Since you are apparently dealing with two polygonal windows, there is no
reason to
invoke the do.call() syntax.
Perhaps the following is what you really want. Note that there are
duplicated points
in the coordinates that you provided (which is not allowed) and that the
coordinates
appear to be given in clockwise order (they must be given in
anticlockwise order).
The following code remedies these problems:
D1 <- unique(data.frame(x=rev(x.coords1),y=rev(y.coords1))
W1 <- owin(poly=D1)
D2 <- unique(data.frame(x=rev(x.coords2),y=rev(y.coords2))
W2 <- owin(poly=D2)
W <- union.owin(W1,W2)
plot(W)
Note that "sp" is not needed for any of the operations in your code (or in
mine).
Finally, please note that gpclib is *NOT* used by the latest version of
spatstat.
It has been replaced by the polyclip package which is not encumbered by the
restrictive license by which gpclib is encumbered. If you are not using the
latest version of spatstat I suggest that you install it forthwith.
HTH
cheers,
Rolf Turner
On 11/02/13 04:30, ASANTOS wrote:
> Dear Members,
>
> I try to create a mask for multiple polygons in spatstat with
> union.owin that is much quicker than the "mask" type, but I have some
> problems that involves list objects and I ask if is posible to use
> union.owin with matrix object? In my example:
>
> require(spatstat)
> require(sp)
> require(gpclib)
>
>
> # Two polygons
> coordinates------------------------------------------------------
>
> #Polygon 1
> x.coords1<-c(371299.9, 371266.4, 371205.6, 371111.8,
> 371047.6, 371018.2, 371014.0,
> 371009.3, 370983.1, 370919.7, 370853.6,
> 370785.6, 370748.8, 370711.8,
> 370687.8, 370696.4, 370785.9, 370885.5,
> 371035.8, 371148.1, 371205.2,
> 371231.7, 371236.5, 371240.3, 371285.8,
> 371326.5, 371397.2, 371417.1,
> 371432.9, 371445.0, 371455.7, 371466.4,
> 371476.6, 371502.6, 371536.0,
> 371550.0, 371546.8, 371528.3, 371470.0,
> 371393.3, 371299.9, 371299.9)
>
> y.coords1<-c(8246589, 8246560, 8246508, 8246428, 8246373,
> 8246349, 8246348,
> 8246352, 8246385, 8246465, 8246551, 8246638,
> 8246685, 8246732,
> 8246764, 8246771, 8246846, 8246932, 8247062,
> 8247160, 8247209,
> 8247230, 8247224, 8247221, 8247160, 8247107,
> 8247016, 8246991,
> 8246967, 8246939, 8246914, 8246892, 8246875,
> 8246846, 8246821,
> 8246809, 8246802, 8246785, 8246735, 8246669,
> 8246589, 8246589)
>
> #Polygon 2
>
> x.coords2<-c(368382.9, 368399.4, 368394.1, 368464.7,
> 368652.2, 368683.7, 368699.7,
> 368714.9, 368714.9, 368714.9, 368860.9,
> 368932.1, 368994.8, 368994.8,
> 369015.7, 369014.5, 369014.5, 369026.8,
> 369045.2, 369020.0, 368952.9,
> 368951.4, 368952.9, 368951.4, 368936.8,
> 368923.4, 368804.1, 368676.0,
> 368178.1, 368182.4, 368202.2, 368216.2,
> 368233.4, 368251.3, 368270.3,
> 368281.7, 368300.6, 368323.8, 368344.8,
> 368364.6, 368374.2, 368382.9)
>
>
>
> y.coords2<-c(8249120, 8249150, 8249153, 8249244, 8249173,
> 8249167, 8249154,
> 8249145, 8249145, 8249145, 8249048,
> 8249026, 8249003, 8249003,
> 8248993, 8248991, 8248991, 8248984, 8248970,
> 8248889, 8248833,
> 8248829, 8248833, 8248829, 8248790, 8248765,
> 8248626, 8248526,
> 8248774, 8248780, 8248815, 8248843, 8248871,
> 8248903, 8248934,
> 8248952, 8248978, 8249013, 8249049, 8249086,
> 8249103, 8249120)
>
>
> # Mask creation
> ---------------------------------------------------------------
>
> mult.poly<-cbind(c(x.coords1,x.coords2),c(y.coords1,y.coords2))
>
> w <- do.call(union.owin,mult.poly)
>
> plot(w)
>
More information about the R-sig-Geo
mailing list