0

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:

all points are stuck at the Northern edge of the leaflet map!

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()

correct projection in basic ggplot2

6
  • It seems inconsistent CRS between points the tiles (map) and what leaflet think expect. (and to me, there are 2 such inconsistencies, so 3 different interpretation of coordinates) Commented Jul 3 at 13:51
  • 2
    In your first line, try using st_transform to project stops into the crs you want - at the moment I think st_set_crs might be just setting it without actually reprojecting the coordinates to match. Commented Jul 3 at 14:26
  • @AndrewGustar I've spent hours trying several combinations of stops <- st_transform(stops, crs = st_crs(3857)) and stops <- st_set_crs(stops, 3857) and st_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?
    – fluegelrad
    Commented 14 hours ago
  • @GiacomoCatenazzi I narrowed down my example to one set of points, so there can only be one mismatch. I set CRS = 3857 for the points and also included the CRS into the leaflet call: leaflet(options = leafletOptions(crs = leafletCRS(crsClass = "L.CRS.EPSG3857"))) Would you know of any other option to specify CRS for the tiles?
    – fluegelrad
    Commented 14 hours ago
  • Tiles: it depends on the source (URL). The tiles names (and location) depends on crs. It is not common not to have a different (so non "WebMecator"). But some publisher may prefer WGS84 Commented 11 hours ago

0

Browse other questions tagged or ask your own question.