Tessellation will not display quadrat totals when plotted (spatstat, R)

Question

I have tessellated an image, and am trying to plot the number of quadrats that appear in each tessellated region.

I have two problems here:

1. My tessellated image will not plot with my gw_ppp points, even with an add = TRUE argument.
2. I am trying to plot the number of quadrants that appear inside each tessellated region, but the image is being plotted with zeros. Also, my gw_ppp points are not being plotted along with the image.

All of the files that can be used to reproduce the error can be found and downloaded in this Google Drive folder (Size: 330Mb).

Here is what I have tried:

library(pacman)
p_load(spatstat,
       dplyr,
       maptools,
       raster,
       sf,
       sp,
       ggplot2)

#Set workplace directory to wherever you downloaded the files from:
setwd(C:\\Users\\Documents)

#Load all data, found in the link above
gw <- read.csv("RileyCNTYGWwells.csv", stringsAsFactors = FALSE)
elev <- raster('elevation.tif')
KS_counties <- st_read("KS_counties.shp")

#Select desired columns from .csv file
gw_sp <- gw %>%
  dplyr::select(LONGITUDE, LATITUDE, WELL_USE, WELL_DEPTH, EST_YIELD) %>% na.omit(gw_sp)

#Convert to spatial dataframe
gw_cor <- st_as_sf(gw_sp, coords=c("LONGITUDE","LATITUDE"),
                crs = st_crs(4326))

#Remove duplicated rows with dplyr's `distinct` function
gw_sp <- gw_cor %>% 
  distinct()

#Omit points outside Riley Co.
riley <- KS_counties %>% 
  filter(name == "Riley")

riley <- st_transform(riley, 4326)

gw_final <- st_intersection(riley, gw_sp)

#Project to a projected CRS, as spatstat is not happy with WGS84
utm14 <- '+proj=utm +zone=14 +ellps=GRS80 +to_meter=0.3048006096012192 +no_defs'

g <- st_transform(gw_final, crs = utm14)

#Filter unwanted columns
g <- g %>% 
  dplyr::select(WELL_USE:EST_YIELD)

#Finally, convert to ppp (RDS file)
gw_ppp <- as(g, "Spatial")
gw_ppp <- as(gw_ppp, "ppp")

From here, we can generate the quadrants:

#Generate the quadrat
q_well <- quadratcount(gw_ppp, nx = 10, ny = 10)

#Plot the quadrats and points
plot(gw_ppp, main = "Riley County Quadrat Well Count", cex = 0.5, pch = "+", cols = "red", legend = FALSE, use.marks = FALSE)
plot(q_well, add=TRUE, textargs = list(cex = 0.8))

These are the quadrant totals that need to show up on the graph (problem 2), but they are not being plotted.

Now, we will read in the elevation raster, reclassify it into 4 categories based on quantiles, and then plot the tessellation together with the gw_ppp points:

#Read in raster and mask to desired county shapefile
elev <- raster("elevation.tif")
riley <- st_transform(riley, crs(elev))
crop_riley <- crop(elev, riley)
mask_riley <- mask(crop_riley, mask = riley)
plot(mask_riley, main = "Riley County Elevation Map")

#Reclassify raster into quantiles
quantile(mask_riley)
elev_zones <- reclassify(mask_riley,
                         c(0, 361.6296, 1,
                           361.6296, 387.6583, 2,
                           387.6583, 403.2133, 3,
                           403.2133, 466.1521, 4))

elev_zones <- ratify(elev_zones)
plot(elev_zones, main = "Elevation Zones")

#Convert to a Spatstat-compatible object
elev_zones <- as.im.RasterLayer(elev_zones)

#Tesselate the image
tes <- tess(image = elev_zones)
plot(tes, main = "Tesselated Elevation Zones")
plot(gw_ppp, add=T, main = "Riley County Quadrat Well Count", cex = 1, pch = "+", cols = "black", legend = FALSE, use.marks = FALSE)

Problem 1 appears after the last line of code is run above. There are no gw_ppp points plotted.

Now, I'm trying to generate the number of quadrants that appear in each tessellated region:

q_elev <- quadratcount.ppp(gw_ppp, tess = tes)
plot(q_elev, main = "Riley County Quadrat Well Count")
plot(gw_ppp, add=T, cex = 1, pch = "+", cols = "black", legend = FALSE, use.marks = FALSE)

Result:

Here is problem 2. These values shouldn't be zeros, and the gw_ppp points are not showing. How can I fix these issues?

Answer 1

tes and everything else derived from the raster is still in lat-long coordinates. Your ppp object is in projected coordinates. plot does not do any transformation to match data with different coordinate systems, so your data do not appear together.

You need to project things to the same coordinate system to display using plot. Which things you project to which CRS really depends on what you plan to do next.