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

[Allar Haav]

Dear Rolf,

Thank you for insightful answer. Following your explanation I found out 
that the key problem was that I didn't specify the simulation in Linhom 
envelope. For some reason I thought that lambda would be overriding the 
envelope creating point process. So, your example 4 was really helpful 
for that matter.

However in order to produce a proper-looking graph with an envelope, I 
had to use reciplambda instead of the regular one. I figured this out 
after seeing that while envelope was calculated and displayed correctly, 
the observed Linhom(r) line did not appear (envelope's "obs" values are 
either 0's or NaN's). Using reciplambda did the trick though, but the 
values on the y-axis are 1e+07 times smaller than they should.  Creating 
random patterns using rpoispp() does not indicate at all that the 
intensity surface should be used reciprocally: the number of points and 
their spatial distribution make sense. Also if using randomly created 
points there is no need for reciplambda, although it works both ways 
showing slightly different results. Another worrying issue, but I'll try 
to look into it.


Thank you and all the best,
Allar


On 23/07/2013 11:07, Rolf Turner wrote:
>
> 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.