[R-sig-Geo] Multiple circular holes with owin function in spatastat

[Rolf Turner]

On 25/07/12 04:13, ASANTOS wrote:
> Dear sig-geoers,
>
>       It's possible to making a multiple circular holes with owin() 
> function in spatstat package if the circles are created in a loop 
> function, for example:
>
> ##
> ##Circle x coordinates
> x<-c(524237,524237,524237,524237,524237,524237,524237,524237,524242
> ,524242,524242,524247,524247,524247,524247,524247,524247,524247
> ,524247,524257,524257,524257,524257,524257,524257,524257)
> ##Circle x coordinates
     You mean "##Circle y coordinates" here.
> y<-c(7978856,7978866,7978876,7978896,7978916,7978936,7978946,7978956,7978956,7978946 
>
> ,7978936,7978956,7978946,7978936,7978916,7978896,7978876,7978866
> ,7978856,7978856,7978876,7978896,7978906,7978916,7978936,7978956)
> ##
>
> CIR<-NULL
> points(a1.past.syn[,1],a1.past.syn[,2])

     What on earth is the forgoing line all about?  It adds nothing to
     the issue of concern and makes your code non-reproducible.
     We don't have "a1.past.syn".
> for (i in 1:length(x)){
>   circle<-disc(radius=5,centre = c(x[i],y[i]))
>   CIR<-rbind(CIR, c(circle))

     Using rbind() here makes no sense at all; actually c(circle) makes no
     sense at all either.
> }
>
> Because I don't know, if I can use CIR's objet in the owin(poly = ... 
> ) function, that uses list() for create this kind of object.

If your discs did not overlap you could just use owin(poly = ...) in the 
following manner:

     CIR <- vector("list",length(x))
     for(i in 1:length(x)) CIR[[i]] <- 
disc(radius=5,centre=c(x[i],y[i]))[["bdry"]][[1]]
     W <- owin(poly=CIR)

***HOWEVER*** your discs *do* overlap, so the foregoing results in an 
error and produces
no window.  So what you have to/can do is:

     CIR <- vector("list",length(x))
     for(i in 1:length(x)) CIR[[i]] <- disc(radius=5,centre=c(x[i],y[i]))
     W <- do.call(union.owin,CIR)

This takes quite a while --- order of half a minute or so on my elderly 
laptop --- and
creates a window of type "mask".

If you really want a window with polygonal boundaries you need to

     * be operating under conditions compatible with the gpclib license

     * have the gpclib package installed

     * enable gpclib by executing

             spatstat.options(gpclib=TRUE)

before executing "W <- do.call(union.owin,CIR)".  It's much quicker than 
the "mask" type
implementation, but the fly in the ointment is that cotton pickin' license!

HTH

     cheers,

         Rolf Turner