[R-sig-Geo] how to limit an image or a funxy to an irregular polygonal window, instead of the whole enclosing rectangle?

Rolf Turner r.turner at auckland.ac.nz
Fri Dec 15 22:20:56 CET 2017


I am probably misunderstanding something (as usual) but I cannot fathom 
what you *expect* the values of the funxy object or the image to *be*, 
outside of the window.

See inline below.

On 16/12/17 08:37, Christopher W. Ryan wrote:
> Using R 3.3.3 and spatstat 1.50-0 on Windows 7.  MWE below.

Thank you for providing a clear and simple example.

> 
> I have a SpatialPolygonsDataFrame of census tract poverty levels in 3
> contiguous counties in the US, called sremsPoverty. I want to use this
> as a predictor in a ppm model. The window for the point pattern is the
> three counties--so an irregular polygon--called sremsWindow.
> 
> I understand ppm predictors need to be an image, a tesselation, a funxy,
> a window, or a single number. So I proceed as follows:
> 
> ### Poverty
> p <- slot(sremsPoverty, "polygons")
> v <- lapply(p, function(z) { SpatialPolygons(list(z)) })
> sat <- tess(tiles = lapply(v, as.owin) )
> pov.f <- as.function.tess(sat, values = sremsPoverty at data$propPoverty)
> 
> Thus pov.f is a spatial function (funxy) that I can use in a call to ppm():
> 
> m1 <- ppm(unique.cases.ppp ~ pov.f)
> 
> pov.f looks as expected when I plot it. But examing the structure of
> as.im(pov.f) it appears it is a 128 x 128 pixel array, with the value of
> the function at all pixels outside of the irregular polygonal window,
> but inside the bounding rectangle, set to NA.

What else could/should they be set to?

> I think this is the cause
> of the NA values I am seeing among the residuals from m1,

I don't think this is the case.  Perhaps more detail is needed; perhaps 
Adrian or Ege will be able to provide more insight.

> and those NA
> residuals are causing me some trouble with model diagnostics such as
> rhohat().

Again I, at least, would need more detail before being able to provide 
any constructive comment.

> How do I constrain the funxy (or the image I can derive from it) to the
> irregular polygonal window, so as to eliminate the NA values outside the
> window but inside the bounding rectangle? Or can I constrain the
> modeling activity of ppm() to the window?

The "modelling activity of ppm()" is *ALWAYS* constrained to the window 
of the pattern to which ppm() is applied.  This is fundamental to the 
the way that ppm() works, and to what a window *means*.

> 
> Thanks.
> 
> --Chris Ryan
> Broome County Health Department
> Binghamton University
> SUNY Upstate Medical University
> Binghamt, NY, US
> 
> MWE:
> 
> x <- c(0, 2.6, 3, 1, 0)
> y <- c(1,2,3,2,1)
> test.window <- owin(poly=list(x=x,y=y))
> plot(test.window)  ## looks as expected
> ## make spatial function
> test.f <- function(x,y){x+y}
> test.funxy <- funxy(test.f, W = test.window)  ## I *thought* this would
> restrict the funxy to the window, but I don't think it does
> plot(test.funxy)  ## looks as expected
> ## but the image from test.funxy is square, with NA values outside the
> polygonal window, which is not square
> test.im <- as.im(test.funxy)
> str(test.im)

Again, what could the values of the image, outside of test.window, 
possibly *be*, other than NA?

> ## my incorrect attempts to restrict the image to the window yields a
> numeric vector
> str(test.im[test.window])

You need to do

      new.im <- test.im[test.window,drop=FALSE]

to get an image rather than a numeric vector.  However the "new.im" that 
you obtain will be identical to test.im:

     > all.equal(new.im,test.im)
     [1] TRUE

> ## or an error message
> window(test.im) <- test.window

You need a capital "W" on Window(test.im) (to avoid an error message). 
But again this won't change anything.

Finally:

set.seed(42)
X <- runifpoint(20,win=test.window)
xxx <- ppm(X ~ test.im)
plot(residuals(xxx))
# No sign of any missing values.

cheers,

Rolf Turner

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



More information about the R-sig-Geo mailing list