Create a map of spatial clusters LISA in R

I would like to create a map showing local spatial cluster of a phenomenon, preferably using Local Moran (LISA).

In the reproducible example below, I calculate the local moran's index using spdep but I would like to know if there is as simple way to map the clustes, prefebly using ggplot2. Help ?

library(UScensus2000tract)
library(ggplot2)
library(spdep)

# load data
data("oregon.tract")

# plot Census Tract map
plot(oregon.tract)

# create  Queens contiguity matrix
spatmatrix <- poly2nb(oregon.tract)

#calculate the local moran of the distribution of black population
lmoran <- localmoran(oregon.tract@data$black, nb2listw(spatmatrix))

Now to make this example more similar to my real dataset, I have some NA values in my shape file, which represent holes in the polygon, so these areas shouldn't be used in the calculation.

oregon.tract@data$black[3:5] <- NA

This is obviously very late, but I came across the post while working on something similar. This uses the rgeoda package, which wasn't out when the question was posted, but it's developed by the GeoDa folks to port some of the functionality of that software to R. sf has also really taken off in the meantime, which makes manipulating spatial data very easy; the rgeoda functions generally expect sf objects.

Like another poster, I'm using the white population instead of Black because the clusters show up better. I converted the original data, with those few observations missing, to sf. rgeoda::local_moran doesn't work when there's missing data, but if you make a copy with the missing observations removed, you can run the analysis and join them back together by ID. Use a right join so you retain all the IDs from the original data, including missing values.

Because this mimics GeoDa, the same colors and labels are stored in the LISA object that local_moran returns. Extract those and use them as the color palette. Because the palette is named, and those names don't include "NA", you can add an NA value to the palette vector, or manually specify a color for NA values to make sure those shapes still get drawn. I made it green just so it would be visible (top left corner).

library(UScensus2000tract)
library(ggplot2)
library(dplyr)
library(sf)
library(rgeoda)

# load data
data("oregon.tract")
oregon.tract@data$white[3:5] <- NA
ore_sf <- st_as_sf(oregon.tract) %>%
  tibble::rownames_to_column("id")

to_clust <- ore_sf %>%
  filter(!is.na(white))
queen_wts <- queen_weights(to_clust)
moran <- local_moran(queen_wts, st_drop_geometry(to_clust["white"]))
moran_lbls <- lisa_labels(moran)
moran_colors <- setNames(lisa_colors(moran), moran_lbls)

ore_clustered <- to_clust %>%
  st_drop_geometry() %>%
  select(id) %>%
  mutate(cluster_num = lisa_clusters(moran) + 1, # add 1 bc clusters are zero-indexed
         cluster = factor(moran_lbls[cluster_num], levels = moran_lbls)) %>%
  right_join(ore_sf, by = "id") %>%
  st_as_sf()

ggplot(ore_clustered, aes(fill = cluster)) +
  geom_sf(color = "white", size = 0) +
  scale_fill_manual(values = moran_colors, na.value = "green") +
  theme_dark()