rspatstat

Kernel density estimation of polygonal window returning NA values

I am using the spatstat package to compute kernel density estimates around an unmarked planar point pattern using a polygonal observation window. The code in the following reprex seems to work as expected up until the density() function which returns an im object with all values of v set to NA whereas I expect the values of v to be the estimated density values.

Any pointers on what I'm getting wrong would be much appreciated!

# Load packages
library(spatstat)
library(sf)

# Create data frame of x and y coords in the Canada Lambert Conformal Conic proj
data <- data.frame(
  x = c(
    8012168.7, -4692979.2, -668427.6, 24334538.4, -12581641.1, 
    2450302.7, -3581952.7, -6652555.7, -9684430.0, 29418984.1
  ),
  y = c(
    1795895, -5823340, 698613, -33489653, 2037623, 
    1182888, -187895, 7096928, -7956464, -26582016
  )
)

# Create a PROJ.4 string for the Canada Lambert Conformal Conic projection
# https://spatialreference.org/ref/esri/canada-lambert-conformal-conic/
clcc_proj <- paste(
  "+proj=lcc",
  "+lat_1=20",
  "+lat_2=60",
  "+lat_0=40",
  "+lon_0=-96",
  "+x_0=0",
  "+y_0=0",
  "+ellps=GRS80",
  "+datum=NAD83",
  "+units=m"
)

# Create a buffer with a radius of 10 km (10000 meters) around each point
data_clcc_buffer <- data |>
  st_as_sf(coords = c("x", "y"), crs = clcc_proj) |>
  st_buffer(dist = 10000)

# Coerce the buffer to an observation window
my_owin <- as.owin(data_clcc_buffer)

# Create a point pattern object
my_ppp <- ppp(data$x, data$y, my_owin)

# Compute kernel density intensity
my_kdi <- density(my_ppp)

# All grid values are NA - why?
all(is.na(my_kdi$v))
Solution

  • Your code creates a window by drawing a circle of radius 10 kilometres around each data point. But the data points are much more widely separated (average nearest-neighbour distance is about 6000 kilometres). At this scale, the window consists of tiny circles around each data point.

    The method density.ppp returns a pixel image with dimensions 128 x 128, by default. A pixel is assigned the value NA if the pixel falls outside the window. The default rule for determining whether a pixel falls inside or outside the window, is to test whether the centre point of the pixel falls inside or outside. In your example, the centres of all the pixels happen to fall outside the window. Therefore the discretised version of the window is empty.

    To avoid this, you could either

    1. specify a larger polygonal window (perhaps by replacing the radius 10 km by a much larger distance)

    2. increase the pixel resolution (using argument dimyx or eps to density.ppp)

    3. first discretise the window using owin2mask

      Window(my_ppp) <- owin2mask(Window(my_ppp), op="cover")

    which will ensure that the polygon is contained inside the pixellated window, and then apply density.ppp.