CODEWITHSUNDEEP
rplotspatialspatstat
Guillem Pocull
Guillem Pocull

Reputation: 25

spatstats: Is there a way to incorporate the GPS accuracy error into the analysis?

I am working on a point-pattern database. I created the spatial analysis using the Linhom function of the Spatstats package on R and it worked perfectly.

However, the points have a GPS accuracy error associated. So I would like to plot this error in the envelopes to compare my observation line, the accuracy error intervals and the envelopes of the complete spatial random pattern.

This is a subset of my dataframe

structure(list(Plant_ID = c("PA0106", "PA0107", "PA0108", "PA0109", PA0110", "PA0111", "PA0112", "PA0113", "PA0114", "PA0115", "PA0116", "PA0117", "PA0118", "PA0119", "PA0120", "PA0121", "PA0122", "PA0123", "PA0124", "PA0125"), error = c(5.7, 3.9, 3, 6.1, 3.6, 6.2, 2.2, 2.3, 2.3, 3.1, 5, 4.2, 6.2, 8.1, 5.7, 3.4, 13.9, 2.4, 3, 2.3), lat = c(4686388.965, 4686405.157, 4686221.337, 4686232.379, 4686373.99, 4686252.22, 4686256.133, 4686373.655, 4686253.222, 4686388.336, 4686387.455, 4686241.46, 4686247.678, 4686387.284, 4686252.61, 4686226.408, 4686254.703, 4686380.872, 4686203.882, 4686215.308), lon = c(423855.312, 423861.0414, 423682.7503, 423725.7385, 423847.7298, 423755.796, 423762.6737, 423847.555, 423773.638, 423854.4161, 423856.1724, 423738.4987, 423747.9986, 423856.7708, 423756.9294, 423717.1068, 423731.609, 423851.6071, 423682.8414, 423711.3172)), row.names = c(NA, 20L), class = "data.frame")

I tried to first create gaussian random simulated points with a number of simulations associated with the GPS accuracy error.

for (i in 1:simulations) {
  
  error_gaussian <- rnorm(nrow(df), 0, df$error/2)
  
  df[, paste0("sim_lat_",i)] <- df$lat + error_gaussian 
  
  error_gaussian <- rnorm(nrow(df), 0, df$error/2)
  
  df[, paste0("sim_lon_",i)] <- df$lon + error_gaussian
  
}

Then I create a random list with the resulting simulations.

random_points_list <- vector(mode = "list", length = simulations)

for (i in 1:simulations) {
  
  x <- df[, paste0("sim_lon_", i)]
  y <- df[, paste0("sim_lat_", i)]
  
  points <- ppp(x, y, window=w)
  random_points_list[[i]] <- rescale(points, 1, unitname = "m")
  
}


random_points_list <- as.solist(random_points_list)

to finally compare it with my real coordinates lat and lon.

real_points <- ppp(df$lon_reproj, df$lat_reproj, window=w)
real_points <- rescale(real_points, 1, unitname = "m")

However, when I plot the function envelope(real_points, Linhom, nsim=simulations, nrank=nrank_alfa, simulate=random_points_list, sigma=bw, eps=1, global=F), it doesn't show the actual behaviour of the simulated points because the envelope interval should get less wide when increasing the radius, isn't it? For example, the first point has an error of 5.7 m, therefore at less than 5.7 m of radius, the envelope should be wider because a lot of simulated points are far away. However, when the radius is 200 m the envelope should be really narrow because the accuracy error is imperceptible.

Here I attach a picture of the resulting plot:

Plot spatstats Lfunction Thank you in advance. I really appreciate your help because I am deeply stuck.

Upvotes: 0

Views: 55

Answers (1)

Adrian Baddeley
Adrian Baddeley

Reputation: 2973

This is an unconventional use of the envelope function but I think it is valid. Usually the envelope of the simulated functions is a significance band for a test of the null hypothesis of (e.g.) complete spatial randomness. However in your example, the simulations are perturbations of the original data attempting to reproduce the GPS uncertainty, so the envelope of the simulated functions is effectively a bootstrap confidence interval for the true value of the inhomogeneous L function. So, the code is fine.

The resulting graph looks OK to me. I think your expectations about the results were not correct. The width of the envelope is related to the amount of variability of the estimate, in absolute terms, not relative to the true value. The L function involves a square root which usually stabilises the variance, so I would have expected the bands to have roughly constant width. It's a bit complicated in this case because you're using Linhom with kernel density estimation of the intensity function, with a fixed bandwidth bw that could be a good choice or a poor choice, so it's difficult to give precise predictions.

Bottom line: results are fine.

Upvotes: 1

Related Questions