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')
EDIT:
ggplot2
plots my points in 3857
without further ado:
stops <- st_as_sf(ausw_input, coords = c("lon", "lat"), crs = 3857)
ggplot(data = stops) +
geom_sf() +
theme_minimal()
st_transform
to projectstops
into the crs you want - at the moment I thinkst_set_crs
might be just setting it without actually reprojecting the coordinates to match.stops <- st_transform(stops, crs = st_crs(3857))
andstops <- st_set_crs(stops, 3857)
andst_crs(stops) <- 3857
. None of them led to the desired leaflet map.ggplot2
easily plots my points (see my edit above), but I'd like to have a map that I can pan and click, just as I could do with polygons for the same area, Wild guess: There's a broken geometry involved? Any hint to check that?leaflet(options = leafletOptions(crs = leafletCRS(crsClass = "L.CRS.EPSG3857")))
Would you know of any other option to specify CRS for the tiles?