I try to plot transit stops on a leaflet map, indicating local clustering ("LISA"), based on Moran's I. The map is embedded in an R markdown file.

I've successfully done that with polygon data, as inspired by this example.

However, when adapting the code to point, the result looks like this:

I've already tried to adjust the CRS, using EPSG:

• 3857, which seems to be the leaflet standard
• 4326, another standard
• 25832, which applies to the area in question (Northern Germany)

stops is my list of transit stops i.e. an sf object of ~7.5k points I want to analyse autocorrelation of the t2 column. My identifier column is h_id_varchar, I use it to assign the labels.

Any suggestions?

The code:

# set CRS
stops <- st_set_crs(stops, 4326)

# create local weights matrix
  nb_loc <- rgeoda::queen_weights(stops)

# Calculate local Moran's I (t2 is the variable I want to investigate)
moran_loc <- rgeoda::local_moran(nb_loc, stops["t2"])

# define labels and colors by the LISA cluster categories
moran_loc_lbls <- lisa_labels(moran_loc)
moran_loc_colors <- setNames(lisa_colors(moran_loc), moran_loc_lbls)

# Define clusters
moran_clustered <- stops %>%
  mutate(cluster_num = lisa_clusters(moran_loc) + 1,  # add 1 bc clusters are zero-indexed
         cluster = factor(moran_loc_lbls[cluster_num], levels = moran_loc_lbls),
         h_id_varchar = as.character(h_id_varchar))  # Ensure that h_id_varchar is character type

# define color palette
pal <- colorFactor(moran_loc_colors, domain = moran_clustered$cluster)


# define LISA label for each location, based on numeric cluster field
moran_clustered <- mutate(moran_clustered, cluster_label = case_when(
  cluster_num == 1 ~ "No significant clustering",
  cluster_num == 2 ~ "<strong>High t2</strong> in this area, <strong>high t2</strong> in the   neighbouring areas",
  cluster_num == 3 ~ "<strong>Low t2</strong> in this area, <strong>low t2</strong> in the neighbouring areas",
  cluster_num == 4 ~ "<strong>Low t2</strong> in this area, <strong>high t2</strong> in the neighbouring areas",
  cluster_num == 5 ~ "<strong>High t2</strong> in this area, <strong>low t2</strong> in the neighbouring areas",
  cluster_num == 7 ~ "<strong>Isolated:</strong>This area has no neighbours",
  TRUE ~ as.character(cluster_num)  # Default label for other cases

))

# assign labels
moran_clustered$label <- paste0(
  "<strong>",
  moran_clustered$h_id_varchar, "<br/>",
  moran_clustered$h_name, "</strong><br/>",
  moran_clustered$cluster_label
) %>%
  lapply(htmltools::HTML)


# render map
leaflet(moran_clustered) %>%
  addProviderTiles('CartoDB.DarkMatter') %>%
  clearShapes() %>%
  addCircleMarkers(
    stroke = TRUE,
    fillColor = ~pal(cluster),
    color = "grey30",
    weight = 1,
    fillOpacity = 1,
    radius = 6,  # Adjust the radius size as needed
    label = ~label,
    layerId = ~h_id_varchar
  ) %>%
  addLegend(position = 'bottomright', pal = pal,
            values = moran_clustered$cluster, title = 'LISA Clusters')