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

[Adrian Baddeley]

Christopher W. Ryan writes:


> 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.

> I think this is the cause of the NA values I am seeing among the residuals from m1,
> and those NA residuals are causing me some trouble with model diagnostics such as
> rhohat().
> 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?

When you convert the tessellation 'sat' into the function 'pov.f' using as.function.tess, the domain of the function is the union of all the tiles in the tessellation 'sat'. After all, this function pov.f(x) works by figuring out which tile of the tessellation the given location x falls inside, and returns the poverty value associated with that tile. So this function is restricted to the union of the tiles given.

When you fit a model to the point pattern 'unique.cases.ppp', the fitting procedure has to place sample points all over the window associated with the point pattern, and evaluate the covariate 'pov.f' at all these sample points. So this window should not be larger than the window where the poverty values are defined - if it is, then you'll get NA's.

I suggest you do something like
        X <- unique.cases.ppp
        W <- intersect.owin(Window(X), Window(pov.f))
        Y <- X[W]
        m1 <- ppm(Y ~ pov.f)

Adrian Baddeley





	[[alternative HTML version deleted]]