Handling missing values in grid-based data for Geary's C

Question

I would like to compute Geary's C or Moran's I indices of spatial autocorrelation for a set of grid-based datasets. I have a bunch of 12x3 grids such as this one:

## Define a 12x3 grid:
loc <- structure(c(0.653, 1.361, 1.392, 0.991, 1.228, 3.882, 1.46, 1.056, 
1.636, 1.535, 1.578, 2.106, NA, 0.801, NA, 1.51, 1.514, 1.846, 
0.791, 1.567, 0.902, 2.011, 1.683, 1.541, 3.437, 1.095, 1.469, 
1.897, 1.866, 1.498, 1.156, 1.342, 1.809, 1.56, 1.341, 1.005), dim = c(12L, 
3L), dimnames = list(c("1", "2", "3", "4", "5", "6", "7", "8", 
"9", "10", "11", "12"), c("1", "2", "3")))

print(loc)

       1     2     3
1  0.653    NA 3.437
2  1.361 0.801 1.095
3  1.392    NA 1.469
4  0.991 1.510 1.897
5  1.228 1.514 1.866
6  3.882 1.846 1.498
7  1.460 0.791 1.156
8  1.056 1.567 1.342
9  1.636 0.902 1.809
10 1.535 2.011 1.560
11 1.578 1.683 1.341
12 2.106 1.541 1.005

As you can see, there may be several missing values in these grids. In R, the package spdep can be used to compute such indices, but it does not handle missing values:

library(spdep)
adj <- cell2nb(nrow = nrow(loc), ncol = ncol(loc),
               type = "queen")
ww <- nb2listw(adj, style = "B")
geary(x = as.numeric(t(loc)),
      listw = ww,
      n = length(as.numeric(loc)),
      n1 = length(as.numeric(loc)) - 1,
      S0 = Szero(ww))

The result we get for that is simply NA.

Theoretically, computing the Geary's C on the grid above is perfectly possible, but it seems that this possibility is not implemented in R. How would be the fastest way to compute an index of autocorrelation for such a grid?

Thanks!