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

[Rolf Turner]

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.