[R-sig-Geo] Question about spatstat: bdist.tiles()

[Rolf Turner]

As Erika  correctly conjectured, the problem boils down to one of rounding
or numerical inaccuracy.  Some of the boundary points of the tiles are 
deemed
to lie outside of the (circular) window in which the pattern is found, 
and these
rejects are not used in calculating the minimum distance from the tile 
to the
window boundary.

A (possible) fix for the problem is to set "check=FALSE" in the call to 
as.ppp()
in bdist.tiles().  This setting results in the points not being 
rejected.  I have
attached the emended code in a file  "bdist.tiles.R".

Place this file in your working directory and do source("bdist.tiles.R") 
before
calling bdist.tiles() --- the new version will mask the one currently in 
spatstat
and thus get called upon, yielding correct answers.

Note that you still don't get DTX[11] being *exactly* 0 (in my toy example).
You get rather 1.82146e-17, and DTX[11] == 0 yields FALSE.  To check
for boundary tiles you need to do something like:

     sapply(DTX,function(x){isTRUE(all.equal(x,0))})

I'm not absolutely certain that my "bug fix" is completely safe, but I'm 
pretty
certain.  I will discuss it with Adrian and if he agrees then this fix 
will be included
in the next release of spatstat.  (And if he disagrees then a different 
fix will be
included!)

     cheers,

         Rolf Turner

On 20/09/12 04:20, Mudrak, Erika [EEOBS] wrote:
> Hello geo list,
>
> I have a question about the bdist function in spatstat when applied to tiles.    The documentation for bdist.tiles says:
> ----------------------
> Details
> This function computes, for each tile s[[i]] in the tessellation X, the shortest distance from s[[i]] to the boundary of the window W containing the tessellation.
> Value
> A numeric vector, giving the shortest distance from each tile in the tessellation to the boundary of the window. Entries of the vector correspond to the entries of tiles(X).
> ---------------------
> *** From what point on the tile does this function measure?  ****
>
> I am interesting is doing Dirichlet (Voronoi) tessellations  (using dirichlet()), but I want to keep separate tiles that hit the window boundary and points that do not.     I thought I would use the bdist.tiles() function and keep only points where bdist>0.
>
> I realize there can be issues with rounding distances, but I came across a situation where there seems to be a bug or an ambiguity for this function.  I have posted an example jpeg at http://mudrak.public.iastate.edu/bdist.tiles_problem.jpg that shows my tessellation, with bdist.tiles values plotted at the point location, and non-boundary tiles (bdist>0) printed in blue. There is a tile on the right that should be a boundary tile but has a bdist of 0.0269, which is larger than many of the other interior tiles.
>
> ** Is there a fix for this bug, or another way to identify boundary tiles? **
>
> Thank you for any suggestions.
>

-------------- next part --------------
bdist.tiles <- function(X) {
  if(!is.tess(X))
    stop("X must be a tessellation")
  W <- as.owin(X)
  switch(X$type,
         rect=,
         tiled={
           tt <- tiles(X)
           if(is.convex(W)) {
             # distance is minimised at a tile vertex
             vdist <- function(x,w) {
               z <- as.ppp(vertices(x), W=w, check=FALSE)
               min(bdist.points(z))
             }
             d <- sapply(tt, vdist, w=W)
           } else {
             # coerce everything to polygons
             W  <- as.polygonal(W)
             tt <- lapply(tt, as.polygonal)
             # compute min dist from tile edges to window edges
             edist <- function(x,b) {
               xd <- crossdist(as.psp(x), b, type="separation")
               min(xd)
             }
             d <- sapply(tt, edist, b=as.psp(W))
           }
         },
         image={
           Xim <- X$image
           # compute boundary distance for each pixel
           bd <- bdist.pixels(as.owin(Xim), style="image")
           bd <- bd[W, drop=FALSE]
           # split over tiles
           bX <- split(bd, X)
           # compute minimum distance over each level of factor
           d <- sapply(bX, function(z) { summary(z)$min })
         }
         )
    
  return(d)
}