I want to calculate the Getis-Ord* statistic for each point in an R simple feature (SF) point object to detect "hotspots" among land cover classes. To do this, I am using the localG function from the spdep package (link) with a set of neighbor and neighbor weights criteria that I have supplied. However, I am getting an error that I am supplying a different number of observations ('Error in localG: XXXX Different numbers of observations'). I'm having a hard time understanding this error and can't find much about it online, so if anyone has ideas about how to interpret and resolve it, I'd love to hear them.
Here is a reproducible example:
# Required libraries
library(dplyr)
library(tidyverse)
library(spdep)
# Data
df <- data.frame(x = c(-113.5974, -113.5947, -113.5920, -113.5893, -113.5866, -113.5839, -113.5812, -113.5974, -113.5947, -113.5920, -113.5893, -113.5866, -113.5839, -113.5812, -113.5974, -113.5947, -113.5920, -113.5893, -113.5866, -113.5839, -113.5812, -113.5974, -113.5947, -113.5920, -113.5893, -113.5866, -113.5839, -113.5812, -113.5974, -113.5947, -113.5920, -113.5893, -113.5866, -113.5839, -113.5812, -113.5974, -113.5947, -113.5920, -113.5893, -113.5866, -113.5839, -113.5812, -113.5974, -113.5947, -113.5920, -113.5893, -113.5866, -113.5839, -113.5812),
y = c(36.99864, 36.99864, 36.99864, 36.99864, 36.99864, 36.99864, 36.99864, 36.99595, 36.99595, 36.99595, 36.99595, 36.99595, 36.99595, 36.99595, 36.99325, 36.99325, 36.99325, 36.99325, 36.99325, 36.99325, 36.99325, 36.99056, 36.99056, 36.99056, 36.99056, 36.99056, 36.99056, 36.99056, 36.98786, 36.98786, 36.98786, 36.98786, 36.98786, 36.98786, 36.98786, 36.98517, 36.98517, 36.98517, 36.98517, 36.98517, 36.98517, 36.98517, 36.98247, 36.98247, 36.98247, 36.98247, 36.98247, 36.98247, 36.98247),
class = c(2, 2, 2, 2, 2, 1, 0, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 0, 1, 0, 2, 2, 2, 1, 1, 2, 0, 2, 2, 2, 0, 1, 1, 1, 1, 2, 2, 0, 0, 1, 1, 2, 2, 1, 0, 1, 1, 1))
wkt <- "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"
sf <- df %>%
st_as_sf(coords = c("x", "y"),
crs = wkt)
## Define neighbors
nb <- cell2nb(nrow=5, ncol=5, type="queen", torus=T)
class(nb)
## Define neighbor weights
list_of_w <- nb2listwdist(nb, sf, type="idw", style="raw",
alpha = 0, dmax = NULL, longlat = NULL, zero.policy=NULL)
## Get Gi*
locG <- localG(sf$class, list_of_w, zero.policy=F, spChk=T)
I expected this code to return a Gi* statistic value (Z-value) for each data point, but instead, I get this error:
Error in localG(sf$class, list_of_w, zero.policy = F, spChk = T) :
Different numbers of observations
What does this error mean and how can I resolve it?
