[R-sig-Geo] spatstat: rhohat and classified raster

Rolf Turner r.turner at auckland.ac.nz
Mon Jun 25 02:09:27 CEST 2012


If I understand your question correctly, it would appear that you are trying
to apply rhohat() in an inappropriate context.  Your covariate seems to have
discrete values; rhohat() is applicable only to continuous covariates.  
Moreover
the theory on which rhohat() is based requires that the spatial cumulative
distribution function of the covariate be *differentiable*.  A (fairly 
minimal?)
sufficient condition for this is that the covariate be differentiable 
with non-zero
gradient.

See:

     "Nonparametric estimation of the dependence of a spatial point process
      on spatial covariates" by  Adrian Baddeley, Ya-Mei Chang, Yong 
Song and
      Rolf Turner, in Statistics and its Interface, vol. 5, no. 2, 2012, 
pp. 221-236.

     cheers,

         Rolf Turner

On 25/06/12 08:15, Yury Ryabov wrote:
> Hi!
>
> I have an issue with interpretation of the results of "rhohat" function
> (spatstat package). I use a classified raster as a covariate to my point
> process. A raster is a land cover image with 300 m resolution. There are
> several classes with values like 10, 20, 50, 100 etc. - all are multiple of
> 10.
>
> The result of rhohat (http://image.bayimg.com/fapcmaadf.jpg) demonstrates
> two peaks at values 135 and 165. Despite the fact that all raster values
> are multiple of 10 and so there are no such values as 135 and 165, there
> are also no such values as 130, 140, 160 or 170 (see histogram:
> http://image.bayimg.com/fapcpaadf.jpg).
>
> So the question is what did I get and why? And is it correct to use a
> classified raster as a covariate in "rhohat"?
>
> Here is the code (P is a ppp):
>
>> gl <- readAsciiGrid("/some_path/land_cover.asc", as.image=TRUE,
>>       proj4string= CRS("+proj=tmerc +lat_0=0 +lon_0=33 +k=1 +x_0=6500000
>> +y_0=0 +ellps=krass
>                              +towgs84=23.92,-141.27,-80.9,-0,0.35,0.82,-0.12
>> +units=m +no_defs"))
> glim <- as.im(gl)
>> plot(rhohat(P, glim, covname="LandCover"),
>>   xlab= "Global Land Cover",
>>   legendpos = "topleft",
>>   legendargs=list(bg="transparent"),
>>   main = NULL)
>
> Best regards,
> Yury Ryabov
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo



More information about the R-sig-Geo mailing list