How to add shapefiles and raster files as covariates in spatstat kppm function

Question

I encountered various errors while running the spatstat package. Here is a summary of my data: a shapefile of landslide point events, a watershed shapefile as the observation window, and several covariates (raster files (*.tif): elevation, slope, aspect, and etc, and vector files: landcover$class_name and lithology$FormUnitNa). I need to run several "cox and cluster" models to see which model is best for my data and identify which landcover classes and lithology type are statistically valuable in my analysis. Finally, use a prediction function to generate a landslide susceptibility map for the entire watershed.

    # Install required packages (if not already installed)
    install.packages(c("ggplot2", "sf", "terra", "spatstat", "raster", "stars"))
    
    # Load libraries
    library(ggplot2)
    library(sf)
    library(terra)
    library(spatstat)
    library(raster)
    library(stars)
    
    # Read the shapefiles (ppplandslide is a marked point pattern)
    shp_lslide <- st_read("landslide_v1.shp")
    shp_wshed <- st_read("Watershed_v1.shp")
    shp_litho <- st_read("lithology.shp")
    
    plot(as.layered())
    
    summary(shp_lslide)
    summary(shp_wshed)
    
    # Create a simplified map using ggplot
    ggplot() +
      geom_sf(data = shp_wshed, fill = "lightblue", color = "blue", alpha = 0.5) + 
      geom_sf(data = shp_lslide, color = "red", size = 2) +
      labs(title = "Landslides within Watershed Boundary") +
      theme_minimal() +
      theme(legend.position = "none")
    
    # Check geometry types (optional)
    st_geometry_type(shp_lslide, by_geometry = FALSE)
    st_geometry_type(shp_wshed, by_geometry = FALSE)
    
    # Convert the landslide centroids to a 'ppp' (point pattern) object
    # First, convert the sf object to a SpatialPoints object, then to ppp
    lslide <- as.ppp(st_as_sf(shp_lslide))
    
    # Convert the watershed shapefile to an 'owin' (observation window) object
    wshed <- as.owin(st_as_sf(shp_wshed))
    
    # Subset the point pattern to include only points within the watershed boundary
    lslide_within_wshed <- lslide[wshed]
    
    lslide_within_wshed <- unmark(lslide_within_wshed)
    # ---- Adding the Elevation Covariate ----
    
    # Load the raster layer (elevation, aspect)
    elevation <- read_stars("DEM_Resample_v2.tif")
    aspect <- read_stars("Aspect_v2.tif")
    slope <- read_stars("Slope_reclass_v3.tif")
    planCurva <- read_stars("PlanCurvature_v2.tif")
    proCurva <- read_stars("ProfileCurvature_Reclass_v1.tif")
    landCover <- read_stars("Land_Cover_v3.tif")
        
    #Convert raster to pixel image
    elevation <- as.im(elevation)
    aspect <- as.im(aspect)
    slope <- as.im(slope)
    planCurva <- as.im(planCurva)
    proCurva <- as.im(proCurva)
    landCover <- as.im(landCover)
    
    plot(elevation)
    plot(slope)
    plot(landCover)
    
    fit1 = kppm(lslide_within_wshed, ~elevation + aspect - 1 + slope -1 + 
              planCurva -1 + proCurva -1 + landCover -1, "Thomas")
    fit2 = kppm(lslide_within_wshed, ~landCover -1, "Thomas")
        
    summary(fit1)
    summary(fit2)

    fitted_intensity <- predict(fit1)
    plot(fitted_intensity, main = "Landslide susceptibility map")

Answer 1

You haven't said what errors were encountered, so I will have to guess.

Your code seems appropriate up to the first call to the function kppm. I am guessing that you did not encounter any problems converting the shapefile data to pixel images of class "im" using as.im. You could check this by typing summary(elevation), is.im(elevation) and so on.

I am guessing that you encountered problems when calling kppm. The call to kppm is rather strange. Instead of

fit1 <- kppm(lslide_within_wshed, ~elevation + aspect - 1 + slope -1 + 
          planCurva -1 + proCurva -1 + landCover -1, "Thomas")

the usual syntax is

fit1 <- kppm(lslide_within_wshed ~ elevation + aspect + slope + 
          planCurva + proCurva + landCover -1, "Thomas")

There is no need for a comma after the name of the point pattern dataset; the point pattern is the left side of the model formula.

The symbol -1 does not need to be included more than once; this is redundant.

Most importantly, we need to understand the meaning of this model formula. For a detailed explanation about model formulas, see Section 9.3 of the spatstat book. In brief, the right hand side of the formula is an expression for the logarithm of the intensity of the point process. Thus the model formula is equivalent to assuming that the log intensity is of the form

log lambda(x) = a_1 * elevation + a_2 * aspect + ... + a_6 * landCover

where a_1, ..., a_6 are coefficients that will be estimated by kppm. (For the moment, I am assuming that all these variables are numeric, not factors). That is, the model states that

lambda(x) = exp(a_1 * elevation + .... + a_6 * landcover)

Therefore the intensity is an exponential function of each of the predictor variables (assuming they are all numeric).

There is no intercept term (constant) on the right hand side, because the model formula includes the symbol -1. It seems strange to exclude the intercept term (if all the variables are numeric). I guess that landCover may be a factor, and in that case, the -1 would ensure that you get a separate coefficient for each level of the factor, which would make sense.

Various kinds of errors can occur. The calculation can become numerically unstable if the predictor variables take very large or very small values. If one of the predictors is a factor, but some levels of the factor are not observed in the data, then the coefficients for these levels will not be estimable (and you should get NA values). It's hard to say more without knowing what error messages you encountered.

Some tips for trouble-shooting:

• use ppm instead of kppm to identify the problem. These two functions use the same mechanism to estimate the intensity function, and ppm is faster.

• try adding the variables to the model one-at-a-time to see when the errors start to happen:

  fit1 <- ppm(lslide_within_wshed ~ elevation)

  fit2 <- ppm(lslide_within_wshed ~ elevation + aspect)

  fit3 <- ppm(lslide_within_wshed ~ elevation + aspect + slope)

• rescale the values of the image predictors to "reasonable" values (bearing in mind that the calculation involves sums of squares of these values). You could use scaletointerval to do that automatically.

• rescale the spatial domain so that the coordinates are not huge numbers (for example, convert metres to kilometres) using rescale

• if you're still having trouble, find the smallest reproducible example of your problem, and post it to SO.