[R-sig-Geo] Using pixel image as spatstat's intensity surface

[Adrian Baddeley]

Allar Haav writes:

> I am currently trying to use spatstat's inhomogeneous cluster analysis
> methods (mostly Linhom), but ran into a problem. Namely, I have a pixel
> image (type "im") with values ranging from 0 to 1 indicating point
> probabilities. When creating and plotting random points with it, this
> intensity surface works as it should.

Intensities are not probabilities.
The intensity is the expected number of points per unit area (possibly spatially-varying).
If the probability of presence of a point is 'p', in a small pixel of area 'a', then the
intensity in that pixel is approximately p/a. 
Thus, to get the correct interpretation in Kinhom, Linhom etc, you need to rescale 
your image pixel values by dividing by the pixel area. 

As a crosscheck, if 'lambda' is an intensity function (stored as a pixel image) that was obtained by fitting
to the point pattern 'X' by some procedure, then integral.im(lambda) should be approximately equal to npoints(X),
because the integral of the intensity is the expected total number of points.

> Linhom(sites, lambda=model)
> # This resulted in:
> "Error in Kwtsum(dIJ, bI, wIJ, b, w = reciplambda, breaks) : Weights in
> K-function were infinite or NA"

This is an important clue that the point pattern called 'sites' apparently includes some points where 
the intensity image called 'model' has a pixel value of either 0 or NA. To check this, extract the 
relevant pixel values at each data point:

       modvals <- model[sites, drop=FALSE]

and inspect them. The inhomogeneous K and L functions attach a weight to each pair of 
data points x_i, x_j proportional to 1/(lambda(x_i) lambda(x_j)) where lambda() is the intensity
function. A value of 0 is illegal because the weight is then infinite, or more insightfully,
because 0 means that "model" thinks it is impossible for a point to occur at that location.
A value of NA is also rejected because that means we don't have the data to compute the inhomogeneous
K function. 

Don't use 'reciplambda' - that is a red herring and will just get everyone confused. 

>  I can also use it in regular
> envelope functions without issues. However, any attempt to use the
> intensity surface with Linhom function results in empty plot.

In calls to envelope() using Kinhom or Linhom, there are two things to consider. 
One is how you want to generate the simulated point patterns that will contribute to the envelope
(this determines the 'null hypothesis' that you are testing).
In your case, you probably want to do something like 
     envelope(..... simulate=expression(rpoispp(mylambda)))
where 'mylambda' is an image containing the intensity values. This will generate realisations
from the inhomogeneous Poisson process with the specified intensity - assuming that is the
null hypothesis.

The other consideration is how you want the simulated points to be weighted in the calculation
of Kinhom or Linhom. For this you could specify the argument 'lambda' that will be passed
from envelope() to Kinhom(). For example if you want to use the same intensity function 'mylambda'
just type
      envelope(X, ..., simulate=expression(rpoispp(mylambda)), lambda=mylambda)

There are some other caveats about the validity of the resulting envelopes (because 'mylambda'
was presumably computed by fitting to the data point pattern 'X') but I'll leave that until later.

Adrian Baddeley

Prof Adrian Baddeley FAA
University of Western Australia  /and/ CSIRO Mathematics, Informatics & Statistics
Mail: <cet.uwa.edu.au>             Skype: adrian.baddeley