Loading Libraries

library(tidyverse) # who doesn't want to be tidy?
library(gt) # for nice tables
library(sf) # for handling geo data
library(here) # reproducible file paths
library(htmltools) # making html labels
library(htmlwidgets) # saving final map
library(leafpop) # pop-up with images
library(leaflet) # mapping

File Folder Names and Loading Data

Here set-up some variables that I use for the file/ folder structure and I read in the spreadsheet.

# folder names
blog_folder <- "posts/2023-08-14-mapping-tombstone"
photo_folder <- "Map"

Reading in the data. This is data created by the previous portion of the project. There are also a bunch of photos in the Map folder. The specific photo that is associated with a given entry is listed in the photo_list column.

tombstones_geo <- readRDS(here(blog_folder, "tombstones_geo.RDS"))

head(tombstones_geo) %>% gt()

Basic Map

simple_map <- leaflet() %>%
  addTiles() %>%
  addScaleBar() %>%
  addCircleMarkers(data = test,
                   label = ~ (paste(complete_name, cemetery_name, sep = " ")))

simple_map

As you can see, no matter how far you zoom, the points never separate. And for what it is worth, the last entry in the test dataframe is shown, indicating the points are added sequentially.

markerCluster Map

Many maps you see on line use a technique called “spiderify” to separate overlapping points. Spiderifying disperses the points as you zoom in.

There are two ways to spiderify a leaflet map. There is a javascript module for spiderifying in can be found here. There is another plug-in called markerCluster, which both clusters and spiderifies points, depending on the zoom level. There is access to this package through the R leaflet package. To access it, use clusterOptions = markerClusterOptions(). For details about the options that can be passed to markerClusterOptions() see the js plug-in documentation.

markerCluster clusters/ groups the data points into a single dot color coded by the number of points it contains (orange, yellow, green) and also displays the number of points in the dot. Single points (not clusters) are displayed in blue or whatever color you set them in the leaflet options. The grouping is set by the number of pixels at the current zoom level (set with maxClusterRadius and defaults to 80 pixels. As you click on a cluster, it zooms in and breaks it up into smaller clusters. (This functionality can be turned off with the parameter zoomToBoundsOnClick.) At the highest zoom, the points are spiderified. This too can be turned off with the parameter spiderfyOnMaxZoom. Both of these features are turned on by default.

I’m going to go back to the 55 data points set just to illustrate this functionality more clearly. (Click on a cluster. It will zoom and split the point as you keep clicking until you end up with just the points with the same coordinates.)

cluster_map <- leaflet() %>%
  addTiles() %>%
  addScaleBar() %>%
  addCircleMarkers(
    data = test,
    label = ~ (paste(complete_name, cemetery_name, sep = " ")), 
    clusterOptions = markerClusterOptions(),
    )

cluster_map

It does separate the points, but it is very hard to click on correct portion of the dot to get the label or pop-up. I think you need to hover on the part of the circle marker that does not overlap with either the cluster marker or the or the other circle markers. This can be solved by making the circle markers very large, so there is more non-overlapping area, but I find this to be ugly.

Luckily, the option spiderfyDistanceMultiplier lets you tailor the distance. I set it to 2, which does put the points km apart (exact distance depends on the zoom). But since they are connected back to the original location with a line, I think it is clear where the points are actually located. When I had other people test the map, the larger displacement was much easier for them to operate the map.

cluster_map_separated <- leaflet() %>%
  addTiles() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addScaleBar() %>%
  addCircleMarkers(
    data = test,
    label = ~ (paste(complete_name, cemetery_name, sep = " ")), 
    clusterOptions = markerClusterOptions(spiderfyDistanceMultiplier = 2)
  ) 

cluster_map_separated

Using sf_jitter

The other option for handling this problem (and my original approach) is to jitter the points using st_jitter(). I think this approach wasn’t as bad as I originally thought. I was guessing which points should have been originally the same location based on cemetery name, but I think probably many of them were in the same cemetery, but with different coordinates to start with.

I found the documentation for this package unclear. Things didn’t behave as I expected. I’m going to go through my work in detail, in case it helps someone else struggling with sf_jitter().

The function st_jitter() takes two parameters, amount and factor. Most examples I’ve seen just use factor, perhaps because everyone found amount as confusing as I did.

Here are the definitions of the parameters.

amount

numeric; amount of jittering applied; if missing, the amount is set to factor * the bounding box diagonal; units of coordinates.

factor

numeric; fractional amount of jittering to be applied

First, it isn’t that hard to find cases where the default value for amount doesn’t work. This can be really disconcerting when recycling working code to other dataframes, and suddenly things are not as expected.

My dataframe test with 3 points at the same location reports a bounding box, but xmin = xmax and ymin = ymax. (The bounding box can be extracted with st_bbox().)

Applying st_jitter() to this dataframe doesn’t produce any sort of errors (and probably should, since then the diagonal of the bounding box is 0.) Instead, the original coordinates are returned with no error or message.

st_bbox(test)

     xmin      ymin      xmax      ymax 
-88.86500  37.98889 -88.86500  37.98889 

And here’s the jittered dataframe. You can see it didn’t jitter it at all.

test_jittered_no_amt <- test %>%
  st_jitter(factor = 1)

test_jittered_no_amt %>% select(geometry)

Simple feature collection with 3 features and 0 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -88.865 ymin: 37.98889 xmax: -88.865 ymax: 37.98889
Geodetic CRS:  WGS 84
                  geometry
1 POINT (-88.865 37.98889)
2 POINT (-88.865 37.98889)
3 POINT (-88.865 37.98889)

Perhaps it is always better to specify amount, so the code will be more robust. This seems reasonable, but I cannot figure out what units amount is in! I think it should be in degrees. The documentation says it is in units of coordinates. The units of the the sf object can be extracted using st_crs().

st_crs(test_for_jitter_inner, parameters = TRUE)$ud_unit

1 [°]

Then, I’d expect amount = 1 and factor = 1 to produce coordinates jittered by about a degree. (The documentation says “For longlat data, a latitude correction is made such that jittering in East and North directions are identical in distance in the center of the bounding box of x.”, so maybe not exactly 1 degree based on the correction.)

test_jittered <- test %>%
  st_jitter(amount = 1,  factor = 1)

test_jittered %>% select(geometry)

Simple feature collection with 3 features and 0 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -88.99865 ymin: 37.08968 xmax: -88.45622 ymax: 38.86458
Geodetic CRS:  WGS 84
                    geometry
1 POINT (-88.99865 38.56987)
2 POINT (-88.45622 37.08968)
3  POINT (-88.5941 38.86458)

Maybe it is the diagonal that is supposed to be 1? This can be calculated with the Pythagorean theorem.

((-89.63298 - -89.32716)^2 + ( 37.39874 - 38.41781)^2)^0.5

[1] 1.063969

So that it, and not at all what I understood from reading the description. The second point of confusion is that I thought that you could specify both factor and amount. (Clearly! As you can see from the example above.) But if you specify amount, then factor doesn’t get used.

In this case, I’d like the points to be jittered by a few meters or so. It is a bit of math to figure out what that is relative to the bounding box diagonal in degrees or some trial and error. On the plus side, the displacement of the dots is independent of the zoom of the map or any other map parameter since it is generated outside the mapping.

So here is the jittered map.

test_jittered <- test %>%
  st_jitter(amount = 0.0001)

map_sf_jitter_all <- leaflet() %>%
  addTiles() %>%
  addScaleBar() %>%
  addCircleMarkers(
    data = test_jittered,
    label = ~ (paste(complete_name, cemetery_name, sep = " "))
  ) 

map_sf_jitter_all

There isn’t any indication that a dot represents multiple points at the more zoomed out scales. Additionally, the zooming experience is very unpleasant. With the clustered option, the map is zoomed and centered as you click through the clusters. There is also a very clear indication (the number in the cluster) that there are multiple data points in the dot. Here, more data points make a darker circle, but this isn’t crystal clear.

Make a nice label formatted with HTML

I’d like a pop-up with biographical information along with the picture, but there is a lot of missing data, so I need to construct the info box carefully or it will be full of NAs.

tombstones_geo %>% select(DOB, DOD)

Simple feature collection with 194 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -89.2775 ymin: 35.97833 xmax: -73.48583 ymax: 44.75056
Geodetic CRS:  WGS 84
First 10 features:
           DOB          DOD                   geometry
1  23 Sep 1877  24 Oct 1899 POINT (-88.69556 38.10361)
2  1 July 1807   2 Feb 1888  POINT (-75.7125 40.87778)
3  17 Apr 1889  19 Nov 1960 POINT (-88.66833 38.06861)
4   5 Aug 1894  1 June 1960    POINT (-88.94 38.19167)
5         1854         1949 POINT (-89.14694 37.60917)
6         1856         1906 POINT (-89.14694 37.60917)
7  27 Feb 1864   3 Dec 1942  POINT (-89.1475 37.60833)
8  29 Sep 1748 19 June 1817 POINT (-75.65028 40.78278)
9  15 Feb 1746  3 June 1823  POINT (-76.0425 41.22528)
10        <NA>         <NA>  POINT (-89.00722 38.1875)

tombstones_geo <- tombstones_geo %>%
  mutate(dob_label = ifelse(is.na(DOB) == TRUE, "", paste0("Born: ", DOB)))

tombstones_geo <- tombstones_geo %>%
  mutate(dod_label = ifelse(is.na(DOD) == TRUE, "", paste0("Died: ", DOD)))

Playing with the map, I see that two entries have the state in both the city and the state column, so I’m going to fix that.

tombstones_geo %>% filter(Surname == "Brown") %>% select(City, State)

Simple feature collection with 2 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -75.7125 ymin: 40.87778 xmax: -75.7125 ymax: 40.87778
Geodetic CRS:  WGS 84
                   City State                  geometry
1 Egypt, Lehigh Co., PA    PA POINT (-75.7125 40.87778)
2 Egypt, Lehigh Co., PA    PA POINT (-75.7125 40.87778)

tombstones_geo <-
  tombstones_geo %>%
  mutate(City = ifelse(City == "Egypt, Lehigh Co., PA", "Egypt, Lehigh Co.", City)
)

Cleaning the City Field

The formatting of the city column is pretty inconsistent. I’m going to clean it up also.

tombstones_geo <-
  tombstones_geo %>%
  mutate(City = str_replace_all(City, "Co\\.", "County")) %>%
  mutate(City = str_replace_all(City, "Co$", "County"))         

In leaflet, it seems that labels require that the html be generated with htmltools::HTML, while popups understand html tags already and can just be passed something like paste("<strong>", complete_name, "</strong>"). Since I was playing around with what info was displayed in label vs. popup, I just rendered everything with htmltools::HTML so I could switch things around without having to modify the text. There is a nice demonstration of the difference between labels and pop-ups on the Dr.Data.King blog.

tombstones_geo <- tombstones_geo %>%
  mutate(
    boxinfo = paste0(
      "<strong>",
      complete_name ,
      "</strong>",
      "<br/><hr>",
      dob_label,
      " ",
      dod_label,
      "<br/><hr>",
      cemetery_name ,
      " in ",
      City,
      " , ",
      State,
      "<br/>"
    ) %>%
      lapply(htmltools::HTML)
  )

Final Map

I’m adding a specific provider tile rather than the default. The label for the point appears when you hover over it and has name and location. There is a pop-up generated in the addCircleMarkers with the biographical info. When addPopupImages() is called, it just appends the photos to that text. The marker pop-up and the pop-up with the picture need to be linked by group.

I couldn’t find a way to directly add text to the addPopupImages() call.

I don’t like the outline around the markers, so I turned it off with stroke = NA and I made the radius = 10.

Next, I call leafpop to add the popup with the images. Documentation for leafpop can be found at CRAN and on the leafpop website. You can add tables, charts, graphs, and images with leafpop. While leafpop can be called within addCircleMarker via popup = popupImage() it will not embed the image within the map. If you intend to save it, then you need to use the separate call to addPopupImages(). Unlike leaflet, the package leafpop does not use the formula notation and requires a “character vector of file path(s) or web-URL(s) to any sort of image file(s).” This is extracted from my dataframe using image_list <- tombstones_geo$photo_list. (Here I use the suffix _list to mean list in the plain English sense of the word, not a list type object.)

I found that specifying the width and max width of the popup was critical. If it were left to the defaults or “too large” then I just got the broken picture icon.

There is also a really weird issue that I can’t figure out that I want to highlight. Sometimes, when I run my quarto document in R, the leaflet map displays with broken pic icons. BUT, rendering the quarto doc does create working pop-ups with the right pics. Other times, both running and rendering produces a working map with pic popups. And even if the map is broken when I “run” instead of “render”, it saves a perfectly working map. I can’t reliably produce either state, and I don’t have any warnings. I couldn’t find anything about this when I googled, so if your map doesn’t work when you are working in quarto (or probably R markdown) try rendering it or saving it.

image_list <- tombstones_geo$photo_list

final_map <- leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addScaleBar() %>%
  addCircleMarkers(
    data = tombstones_geo,
    label = ~ lapply(
      paste("<strong>", complete_name, "</strong>"),
      htmltools::HTML
    ),
    popup = ~ boxinfo,
    clusterOptions = markerClusterOptions(
      spiderfyDistanceMultiplier = 2,
      maxClusterRadius = 50
    ),
    radius = 10,
    stroke = NA,
    group = "group1"
  ) %>%
  leafpop::addPopupImages(
    image = paste0(here(blog_folder, photo_folder), "/", image_list),
    src = local,
    group = "group1",
    width = 400,
    maxHeight = 300,
    maxWidth = 400
  )

final_map

The map can be saved with saveWidget() from htmlwidgets.

# fancy interactive map to upload
saveWidget(final_map, file = "map_to_upload.html")