Interpolation within window using spatstat leaves out one point

Question

I have a data frame of a mapped forest plot, where all tree stems have X, Y coordinates, diameter at breast height (in cm), and survival (0,1) (found here, named "MFP14_surv_forSO.csv": https://www.dropbox.com/sh/t10b53qcobvxlzg/AACZyASgudtFLiZ79QRIjHH_a?dl=0).

I created a spatial point pattern and then a smoothed kernel (im object) of large stem death as a proxy for canopy "gappiness" (or amount of light penetration). I consider large stems as those with a diameter >9cm, and a '0' survival indicates the stem died.

I'd like to interpolate this measure of "gappiness" onto all tree stems. I already know that there are several points in the pattern that lie outside the window, so I exclude those from the analysis using inside.owin() to subset the data frame.

surv14 <- read.csv("MFP14_surv_forSO.csv")
win14 <- owin(poly=list(x=c(0,250,250,225,225,0),y=c(0,0,50,50,100,100)))  #specifying window extent
surv14 <- surv14[inside.owin(surv14[,1],surv14[,2], win14)==TRUE,]  #removing points outside of window
death <- surv14[!is.na(surv14$diam90) & surv14$diam90>9,]  #subsetting only large stems
death <- death[death$surv==0,]  #subsetting only the large stems that died
death.pp <- as.ppp(death,win14)  #creating point pattern from large stem death
death.fun <- Smoothfun(death.pp,sigma=10,edge=TRUE)  #smoothed kernel of large stem death
im <- as.im(death.fun)  #converting smoothed kernel into im object
ext <- im[surv14[,1:2]]  #yields pixel values of gappiness for each stem
surv14 <- cbind(surv14,ext)  #adding gappiness measure to data frame

But when I interpolate, it is leaving out one of the points - there are 4873 observations in surv14 and only 4872 in the interpolation, ext. So when I try to bind the measure of gappiness to my data frame using cbind I get the following error message:

Error in data.frame(..., check.names = FALSE) : 
arguments imply differing number of rows: 4873, 4872

I don't know how to figure out which point it's leaving out, and why. Any guidance would be greatly appreciated!

Answer 1

The reason why this fails is a bit complicated.

I will try to go through your code and explain what's going on:

library(spatstat)
surv14 <- read.csv("MFP14_surv_forSO.csv")
win14 <- owin(poly=list(x=c(0,250,250,225,225,0),y=c(0,0,50,50,100,100)))  #specifying window extent
surv14 <- surv14[inside.owin(surv14[,1],surv14[,2], win14)==TRUE,]  #removing points outside of window
death <- surv14[!is.na(surv14$diam90) & surv14$diam90>9,]  #subsetting only large stems
death <- death[death$surv==0,]  #subsetting only the large stems that died
death.pp <- as.ppp(death,win14)  #creating point pattern from large stem death

At this point you have a marked point pattern with two different marks: diam90 and surv, and to avoid complications I will only keep diam90 (since we know surv is 0 for all the points):

marks(death.pp) <- marks(death.pp)$diam90

Continuing your code you next create death.fun which is an R function of two arguments x and y and it returns the kernel smoothed tree diameter at the location x,y which must be inside the polygonal window win14 defined above.

death.fun <- Smoothfun(death.pp, sigma=10, edge=TRUE)  #smoothed kernel of large stem death

Now you convert this function to a pixel image with default resolution 128x128 pixel over the enclosing rectangle.

im <- as.im(death.fun)  #converting smoothed kernel into im object

On a simple plot with the original window overlayed the digital approximation looks fine:

plot(im, main = "")
plot(win14, add = TRUE, border = "green", lwd = 3)

However, you cannot perfectly represent the original window in this resolution. Point number 1850 of surv14 is inside the original window but not the digital approximation:

inside.owin(surv14$X[1850], surv14$Y[1850], win14)
#> [1] TRUE
inside.owin(surv14$X[1850], surv14$Y[1850], Window(im))
#> [1] FALSE

Then you are subsetting a pixel image by a two-column data.frame which is interpreted as a list of points and the pixel value for each point inside the window of the image is returned. Hence only 4872 points are returned:

ext <- im[surv14[,1:2]]  #yields pixel values of gappiness for each stem
length(ext)
#> [1] 4872

The simple fix for this is to use death.fun directly and skip the pixel approximation:

ext <- death.fun(surv14$X, surv14$Y)
#> Warning in ppp(X$x, X$y, window = win, check = check): data contain
#> duplicated points
length(ext)
#> [1] 4873

Then your final line of code should work:

surv14 <- cbind(surv14,ext)  #adding gappiness measure to data frame

Please note that I have not considered whether it is a good idea to use this measure of gapiness here. I will leave that to you.

Answer 2

AJ Brown asked:

I don't know how to figure out which point it's leaving out, and why.

In the command ext <- im[surv14[,1:2]] you are using [.im to extract pixel image values at the locations with coordinates surv14[,1:2]. The help for [.im says that it has an argument drop that specifies whether locations that lie outside the pixellated window should be omitted (drop=TRUE) or should return an NA value (drop=FALSE). The default is drop=TRUE so that if you have any locations outside the pixellated window the resulting vector ext will be shorter than expected. So if you do

ext <- im[surv14[,1:2], drop=FALSE]

you can then use is.na(ext) to find the offending points.

You don't need the pixellation in this particular example, so it's better to just use the Smoothfun, as Ege advised.