Using pixel image as spatstat's intensity surface

Dear Allar,

It is impossible for me to tell from the information that you have 
provided what is going
wrong.  Basically it seems to me that the code you used will work "in 
general".  It is strange
that you are getting *empty* plots.   A toy example (see below) 
indicates that "empty plots"
should not result.

I can see why there might be problems with the envelope, but the plot 
shouldn't be
*empty*.  At the very least you should get a plot of Linhom(sites).

Toy example:

W <- disc(0.5,c(0.5,0.5),npoly=6)
A  <- as.im(function(x,y){0.1*cos(2*pi*(x^2+y^2))^2},W=W)
i    <- summary(A)$integral
B   <- eval.im(100*A/i)
set.seed(42)
X  <- rpoispp(B)[W]
E1 <- envelope(X,Linhom,lambda=A,nsim=19)
E2 <- envelope(X,Linhom,lambda=B,nsim=19)
E3 <- envelope(X,Linhom,lambda=A,simulate=expression(rpoispp(A)),nsim=19)
plot(E1)
plot(E2)
plot(E3)

The envelopes E1 and E2 are silly looking, being inappropriate, but 
they are
certainly not empty.

The envelope E3 yields a non-empty plot, although there is no envelope 
as such.
This is due to the fact that the "simulate" expression yields nothing 
but empty
patterns (from which no L function estimate can be calculated).

I suspect that what you actually did was something like the E3 
calculation rather
than what you told us you did.

What you probably should have done is something like:

E4 <- envelope(X,Linhom,lambda=B,simulate=expression(rpoispp(B)),nsim=19)

which produces a sensible envelope.

What Linhom() does for you is to allow you to assess whether "sites" 
arises from an
inhomogeneous Poisson process or if there is evidence of interaction 
in the data.  To
apply Linhom() you must know (or have an idea of) the underlying 
intensity of "sites".
To use an intensity such as "model" which you apparently *know* would 
generate fewer
points than there are in "sites" makes no sense.  Use an intensity 
which you actually
believe to be the underlying intensity for "sites" (or which you at 
least believe to be a
*plausible* candidate for this intensity).

    cheers,

        Rolf Turner

P. S.  The error that you got from applying Linhom() directly to 
"sites" with lambda=model
might have arisen because "model" has some zero values.  I *think* 
this will be a problem
only if it has zero values at one or more points of "sites". Haven't 
experimented with
this.

        R. T.

On 23/07/13 13:23, Allar Haav wrote:

Hi all,

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. 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. I suspected that this might have something to do with the 
integral of the surface, therefore made a crude transformation so 
that it would be approximately the number of points under study (see 
below), but of course it didn't change anything. Generally I still 
feel that the crux lies in the intensity surface's values, but I 
can't quite put my finger on it.

A simplified code I used is as follows:

area <- as.owin(area_polygon)
sites <- as.ppp(as(readShapeSpatial('sites.shp')))
sites$window <- area
model <- as.im(read.asciigrid('model.asc'))
sites.n <- sites$n
integral <- summary(model)$integral

# This creates a nice plot with simulation envelopes
plot(envelope(sites, Lest, nsim=19, 
simulate=expression(rpoispp(eval.im(model*(sites.n/integral ))))))

# These produce only empty plots
plot(envelope(sites, Linhom, nsim=19, lambda=model ))
plot(envelope(sites, Linhom, nsim=19, lambda=eval.im(model 
*(sites.n/integral ))))

# A further check without envelope
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"

That last error message left me a bit perplexed. Indeed, I have NA 
values in the dataset as the study area is not perfectly rectangular. 
But this should be a non-issue anyway.
Frankly I'm running out of ideas. I've been thinking of alternative 
approaches, e.g. creating random points from the intensity surface 
and fitting a model based on it using "ppm", perhaps using some 
covariates. Yet this seems actually an awful idea considering the 
probable loss of accuracy and the fact that basically the model 
already exists.

As a side note, a bit of googling on the issue revealed an older post 
about inhomogeneous process analysis 
(http://r-sig-geo.2731867.n2.nabble.com/memory-allocation-problem-with-envelope-function-in-spatstat-td5029218.html). 
There a command equivalent to my first envelope plot was used: 
Lin=envelope(SPECIES.ppp,Linhom,nsim=5,simulate=expression(rpoispp(lambda)), 
! + correction="border")
That made me think that I might be overthinking the problem and I've 
already found the solution with first try. However I am still 
skeptical as the difference between Lest and Linhom does not lie 
merely in simulated patterns - weights are also to be considered. If 
so, why was this solution used?

Any suggestions or ideas? All sort of feedback would be appreciated.