[R-sig-Geo] silly question on quadratcount in spatstat ?

[Adrian Baddeley]

Alexendre Villers writes: 

> library(spatstat)
> data(humberside)
> plot(humberside)
> qH <- quadratcount(humberside, 5, 5)
> plot(qH, add=TRUE, col="blue", cex=1.5, lwd=2)

> There is obviously something wrong with the display of the count 

What is wrong with the display?

Since the window for 'humberside' is an irregular shape, when you issue the command 'quadratcount(humberside, 5, 5)' the rectangle surrounding this region is first divided into a 5 x 5 array of rectangles, and then these rectangles are intersected with the irregular region to form the quadrats. At this stage there are 25 quadrats, some of which may be empty because the rectangle did not intersect the irregular window. The number of data points in each quadrat is counted (yielding 25 numbers). Also the area of each quadrat is counted (yielding 25 numbers, some of which may be zero).

When the object qH is plotted, the plot method shows only the quadrats that have nonzero area (because they intersected the window). For each of these quadrats, the count is displayed as a numeral at the 'centre' of the quadrat. The centre is computed by finding the approximate incentre of the quadrat (the centre of the largest disc inscribed in the quadrat). 

In the case of 'quadratcount(humberside, 5, 5)' some of the quadrats are very small pieces. In these small quadrats, the plot method (plot.quadratcount) has some difficulty in finding a comfortable place to put the numeral that represents the count. 

> qH gives a contingency table of counts on a 5 * 5 grid but there are 
> only 20 cells lying inside the owin. Why do we get value for cells 
> outside our observation window ?  

Well, this is just a convention, so that quadratcount(X, 5, 5) always results in a contingency table with 25 cells. 

> I'm sure there is a way of getting a vector for the cells 
> observed it but I couldn't figure how. 

The easiest way to get all the information you want is to use quadrat.test(). 
This removes quadrats that do not intersect the window. 
You can extract the observed counts and expected counts as follows:

           qtH <- quadrat.test(humberside, 5, 5)

           observed.counts <- qtH$observed
           expected.counts <- qtH$expected

If you want the actual quadrats, 

           quadratsH <- tiles(as.tess(qtH))

If you want the areas of the quadrats,

           quadrat.areas <- unlist(lapply(quadratsH, area.owin))

regards
Adrian Baddeley