how to limit an image or a funxy to an irregular polygonal window, instead of the whole enclosing rectangle?

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

how to limit an image or a funxy to an irregular polygonal window, instead of the whole enclosing rectangle?

Christopher W. Ryan
Using R 3.3.3 and spatstat 1.50-0 on Windows 7.  MWE below.

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@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?

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)

## my incorrect attempts to restrict the image to the window yields a
numeric vector
str(test.im[test.window])
## or an error message
window(test.im) <- test.window

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: how to limit an image or a funxy to an irregular polygonal window, instead of the whole enclosing rectangle?

Rolf Turner

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

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: how to limit an image or a funxy to an irregular polygonal window, instead of the whole enclosing rectangle?

Christopher W. Ryan
Thanks Rolf. I'm going to have to reflect more on my code and my data,
to understand better what is going on.

Obviously this won't help you much, without having access to all my data
and preceeding code, but the error message that is tripping me up is:

> rhohat(m12, pov.f)
Error:  the fitted intensity contains NA values

And yet,

table(is.na(fitted(m12)))
FALSE
  876

The predicted intensity, however, contains many NA values:

table(is.na(predict(m12)$v))
FALSE  TRUE
10379  6005

I try to force predictions only within my window by specifying locations
(which I think requires a binary mask window) but get the same result:

> table(is.na(predict(m12, locations = as.mask(sremsWindow))$v))
FALSE  TRUE
10379  6005

Does rhohat use fitted values (at quadrature points) or predicted values
(on a 128 x 128 pixel grid within the window)? Top of p. 415 in your
book Spatial Point Patterns: Methodology and Applications wtih R seems
to indicate the latter, while the error message from my rhohat command
above speaks of fitted values.  And how is a rectangular 128 x 128 grid
fit in an irregularly-shaped polygonal window?  Maybe that's how NA
predicted values arise--pixels outside an intra-window rectangular grid
but still inside the window?


And I can see now that no residuals from the model are missing:

> table(is.na(residuals(m12)$val))
FALSE
  876

All the NA's in my predicted values *around* my window, but inside the
bounding rectangle, led me down the garden path.

The origin of most of my predictors, such as pov.f, are shapefiles from
the US Census Bureau, with a discrete value of poverty level for each
census tract. So a tesselation of my window, really.  Through much
wrangling (possibly poorly-done) I was able to turn each predictor into
a funxy--therefore they are essentially step functions, constant across
a census tract and with abrupt changes at census tract boundaries.  I
notice that rhohat calls for a continuous covariate. Could that be an
issue?  Although, I have one predictor that is a continuous distfun, and
I get the same error message from rhohat with that.


Thanks.

--Chris Ryan

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

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: how to limit an image or a funxy to an irregular polygonal window, instead of the whole enclosing rectangle?

Rolf Turner
On 16/12/17 11:43, Christopher W. Ryan wrote:
> Thanks Rolf. I'm going to have to reflect more on my code and my data,
> to understand better what is going on.
>
> Obviously this won't help you much, without having access to all my data
> and preceeding code,

Too true!

> but the error message that is tripping me up is:
>
>> rhohat(m12, pov.f)
> Error:  the fitted intensity contains NA values

No idea what is causing that to happen, and it's impossible to work it
out without having your data.  (Point pattern and the predictor image.)

I can't find that error message anywhere in the spatstat code.  Is it
possible that you have some other rhohat() function (maybe from some
other package) hanging around and getting in the way?  Using traceback()
*might* provide some insight.

>
> And yet,
>
> table(is.na(fitted(m12)))
> FALSE
>    876
>
> The predicted intensity, however, contains many NA values:
>
> table(is.na(predict(m12)$v))
> FALSE  TRUE
> 10379  6005

This is a *COMPLETE RED HERRING* !!!

Of course you will get NAs from this.  The predicted intensity is an
*image* on the window of the point pattern to which the model is fitted.
All pixels within the bounding box of that window, but outside of the
actual window, have pixel value NA.
>
> I try to force predictions only within my window by specifying locations
> (which I think requires a binary mask window) but get the same result:
>
>> table(is.na(predict(m12, locations = as.mask(sremsWindow))$v))
> FALSE  TRUE
> 10379  6005

No, no, no!!!  The predictions are always "within your window".  That's
*why* they are NA at pixels outside that window.

>
> Does rhohat use fitted values (at quadrature points) or predicted values
> (on a 128 x 128 pixel grid within the window)? Top of p. 415 in your
> book Spatial Point Patterns: Methodology and Applications wtih R seems
> to indicate the latter, while the error message from my rhohat command
> above speaks of fitted values.

I don't understand where that error message is coming from, so I don't
get what it is on about, but essentially rhohat() looks at values of the
predictive covariate and of the fitted intensity at all pixel centres
within the window.  The 128 x 128 grid is the default here, but can be
changed (in the call to rhohat()).

> And how is a rectangular 128 x 128 grid
> fit in an irregularly-shaped polygonal window?  Maybe that's how NA
> predicted values arise--pixels outside an intra-window rectangular grid
> but still inside the window?

The rectangular grid is over the bounding box of the window.  (Which is
rectangular!!!) Only pixels whose centres are within the window have
non-NA values.

>
>
> And I can see now that no residuals from the model are missing:
>
>> table(is.na(residuals(m12)$val))
> FALSE
>    876
>
> All the NA's in my predicted values *around* my window, but inside the
> bounding rectangle, led me down the garden path.

Indeed.

>
> The origin of most of my predictors, such as pov.f, are shapefiles from
> the US Census Bureau, with a discrete value of poverty level for each
> census tract. So a tesselation of my window, really.  Through much
> wrangling (possibly poorly-done) I was able to turn each predictor into
> a funxy--therefore they are essentially step functions, constant across
> a census tract and with abrupt changes at census tract boundaries.  I
> notice that rhohat calls for a continuous covariate. Could that be an
> issue?

I don't think so.  The rhohat() function works by smoothing, and
smoothing a step function doesn't really make sense.  But I think that
rhohat() should still give you an answer, even though it's a crappy answer.

> Although, I have one predictor that is a continuous distfun, and
> I get the same error message from rhohat with that.

Don't think that I can come up with any useful suggestions as to how to
proceed from here.  Sorry.

cheers,

Rolf Turner

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

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: how to limit an image or a funxy to an irregular polygonal window, instead of the whole enclosing rectangle?

Adrian Baddeley-3
In reply to this post by Rolf Turner
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@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]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo