[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