[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
Sat Dec 16 06:13:53 CET 2017
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
More information about the R-sig-Geo
mailing list