This week, The DataLab made a map of “amenities” that are near Reed. The data comes from OpenStreetMap (OSM), which has a huge amount of useful and free data, but also some very interesting classifications for things! Read on to find out more about why that is and how we made the map.

OSM is a free, editable map of the world built by volunteers. It’s basically Wikipedia for maps. All of its data comes from user-contributed tags, so anyone can add and label features. These are often things like classifying things as residential or commercial with a tag like “landuse=residential” or describing natural features like “waterway=river”.

But people can tag whatever they like, and a lot of fun things are tagged under “amenity”. There are many practical things like fire stations, but there are also tags for things like “icecream” and, my favorite, “dog_toilet”. Sadly, most of the dog toilets are in Europe and there are no tagged locations of dog toilets near Portland, but because anyone can add tags, you could change that!

Where do you find the data?

All of the data comes from OpenStreetMap and was accessed directly in the R script.

For the print map, I pulled the data live from the Overpass API, which is a free service for querying OSM. The {osmdata} package lets you draw a “bounding box” (just the latitude and longitude boundaries for your map) and then you can query for every feature with a certain tag inside that box. That worked fine for the small area around Reed, but Overpass was unreliable for anything larger than the small box.

# Bounding box vaguely around inner Portland to I-205 on the east
bb <- c(-122.72, 45.46, -122.55, 45.56)

For the interactive map, I wanted the user to be able to zoom out a bit and look at areas around broader Portland, and Overpass kept giving me 502 and 504 errors. This means the public servers were overloaded and big requests were failing. So instead, I downloaded a whole OSM extract of the Portland metro area from BBBike using the {osmextract} package. That file is large, so it isn’t included in the Github repo for this project, but a cached rds file is, so if you’re looking only for the same things I wanted, the repo works well. If you want to expand beyond the amenities I chose to look at, you’ll need to extract the larger file yourself.

# my local file
extract_file <- "/Users/griffinj/Documents/GIS_files/osm_extracts/bbbike_Portland.osm.pbf"

# the cache it creates
cache_file <- "amenities/data/amenities_raw_inner_pdx.rds"

I pulled amenities using the add_osm_feature() function and decided to also add another feature to pull dispensaries, which we’re classified differently as “cannabis shops”. Some of the locations that were pulled were mapped as single points and others were polygons in the shape of the outline of their buildings. I just wanted to plot everything as points, so I converted the polygons to “centroids” (the center point of the polygon) so that everything could be plotted the same way.

# --- Pull amenities from OpenStreetMap -----------------------------------
# Source: OpenStreetMap contributors via the Overpass API (osmdata package).
# Amenities are points, but some are tagged on building footprints, so we
# keep the point and polygon layers and reduce polygons to their centroids.

# Points, plus polygons (buildings) reduced to their centroids.
as_points <- function(osm) {
  polygons <- if (is.null(osm$osm_polygons)) NULL else st_centroid(osm$osm_polygons)
  bind_rows(osm$osm_points, polygons)
}

amenities_osm <- opq(bbox = bb) |>
  add_osm_feature(key = "amenity") |>
  osmdata_sf()

# # Dispensaries are tagged as shops, not amenities, so they need their own query.
dispensaries_osm <- opq(bbox = bb) |>
  add_osm_feature(key = "shop", value = "cannabis") |>
  osmdata_sf()

amenities <- as_points(amenities_osm) |>
  filter(!is.na(amenity)) |>
  transmute(name, amenity)

dispensaries <- as_points(dispensaries_osm) |>
  filter(shop == "cannabis") |>
  transmute(name, amenity = "dispensary")

amenities <- bind_rows(amenities, dispensaries) |>
  st_transform(4326)

How do you make the map?

For the plain print map, I needed a pretty sparse background, so I just added streets from OSM and a shapefile of the outline of Reed.

# add the streets
streets_osm <- opq(bbox = bb) |>
  add_osm_feature(key = "highway", value = c("motorway", "trunk", "primary", "secondary", "tertiary", "residential")) |>
  osmdata_sf()
streets <- streets_osm$osm_lines

# Reed campus boundary 
reed <- st_read("amenities/data/reed_boundary/reed_boundary_WGS.shp", quiet = TRUE) |>
  st_transform(4326)

For the interactive map, I could have a more detailed background and I wanted the user to be able to zoom in and out, so there I used the {leaflet} package to add an ESRI map as the background. {leaflet} is a wonderful package for making interactive maps, so I highly suggest you look into it if you’re going to be working with spatial data. ESRI is an evil overlord that has cornered the market on a lot of GIS software, but they do have a good basic background map. It’s easily added to leaflet plots with just the line addProviderTiles("Esri.WorldGrayCanvas").

The interactive map was actually much easier to make than the smaller black and white map. This is the code that creates the interactive map:

# make the map with leaflet so it's interactive
map <- leaflet(shared, width = "100%") |>
  # provider tiles is what the background of the map is
  addProviderTiles("Esri.WorldGrayCanvas") |>
  # make the Reed boundary distinct
  addPolygons(data = reed, fillColor = "grey50", fillOpacity = 0.35, color = NA) |>
  # add the points
  addCircleMarkers(
    lng = ~lon, lat = ~lat,
    radius = 6, stroke = TRUE, color = "white", weight = 1,
    fillColor = ~pal(type), fillOpacity = 0.9,
    label = ~lapply(paste0(htmlEscape(name), "<br/>", htmlEscape(type)), HTML)
  ) |>
  # center initial map
  setView(lng = reed_center[1], lat = reed_center[2], zoom = 14)

The smaller map took so much more fiddling. By doing these Quest graphics, I’m learning a lot about how hard it is to create useful visualizations when there are a lot of restrictions on style. This map had to be in black and white, so I used shapes to symbolize amenity type, but some places were right next to each other (the places on the corner of Bybee and Milwaukie), so the individual shapes weren’t legible. Making those not overlap was perhaps more trouble than it was worth, but I get stubborn.

To fix the overlap, I grouped any places within 60 meters of each other into a cluster using hclust(). Then I moved each symbol in a cluster onto a small ring around the cluster’s center and drew a thin line back to the true location. That required actually doing some math, so I had to temporarily convert the coordinates to UTM, which is a projection measured in meters instead of degrees. Then I had to actually remember some freaking geometry to get that lines to go to the right places. Here’s the code for the simpler map:

# --- Map -----------------------------------------------------------------
# Keep the most common amenity types and lump the rest into "other".

amenities_plot <- fun_stuff |>
  filter(amenity %in% c("dispensary", "ice_cream", "swingerclub", "pub")) |>
  mutate(type = fct_infreq(amenity))

# Spread out overlapping symbols. Places within `merge_dist` meters of each other
# form a cluster; a cluster's symbols are moved onto a small ring around its center
# and a line points back to each true location.
merge_dist <- 60  # meters
ring_base <- 75   # ring radius = ring_base + ring_step * cluster size, in meters
ring_step <- 15

pts_m <- amenities_plot |> st_transform(26910)  # UTM 10N, meters
xy <- st_coordinates(pts_m)

amenities_plot <- amenities_plot |>
  mutate(
    cluster = cutree(hclust(dist(xy), method = "complete"), h = merge_dist),
    x0 = xy[, 1], y0 = xy[, 2]
  ) |>
  st_drop_geometry() |>
  group_by(cluster) |>
  mutate(
    n = n(),
    angle = 2 * pi * (row_number() - 1) / n,
    radius = if_else(n > 1, ring_base + ring_step * n, 0),
    x1 = mean(x0) + radius * cos(angle),
    y1 = mean(y0) + radius * sin(angle)
  ) |>
  ungroup()

# Back to lon/lat for plotting.
to_lonlat <- function(x, y) {
  st_as_sf(tibble(x, y), coords = c("x", "y"), crs = 26910) |>
    st_transform(4326) |>
    st_coordinates()
}
true_ll <- to_lonlat(amenities_plot$x0, amenities_plot$y0)
shown_ll <- to_lonlat(amenities_plot$x1, amenities_plot$y1)

amenities_plot <- amenities_plot |>
  mutate(
    lon0 = true_ll[, 1], lat0 = true_ll[, 2],
    lon1 = shown_ll[, 1], lat1 = shown_ll[, 2]
  )
moved <- filter(amenities_plot, n > 1)

amenities_map <- ggplot() +
  geom_sf(data = reed, fill = "grey75", color = NA) +
  geom_sf(data = streets, color = "grey80", linewidth = 0.5) +
  geom_segment(data = moved, aes(x = lon0, y = lat0, xend = lon1, yend = lat1), color = "grey30", linewidth = 0.3) +
  geom_point(data = moved, aes(x = lon0, y = lat0), color = "grey30", size = 0.6) +  # true location
  geom_point(data = amenities_plot, aes(x = lon1, y = lat1, shape = type), color = "black", size = 3.5) +
  coord_sf(xlim = bb[c(1, 3)], ylim = bb[c(2, 4)], expand = FALSE) +
  scale_shape_manual(
    name = NULL,
    values = c(dispensary = 16, ice_cream = 17, swingerclub = 18, pub = 15),  # circle, triangle, diamond, square
    labels = c(dispensary = "Dispensary", ice_cream = "Ice cream shop", swingerclub = "Swingers club", pub = "Pub")
  ) +
  theme_void() +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    plot.caption = element_text(color = "grey40"),
    legend.text = element_text(size = 12),
    plot.margin = margin(10, 10, 10, 10)
  )

And that’s for only plotting four amenities! You’ll notice that for doing this map all I needed was good ole ggplot and the geom_sf() command from the package {sf}. That package allows for plotting spatial data within ggplot provided your data objects have either latitude and longitude or a geometry column, if you’re plotting shapes. For a long time, I thought the “sf” in the package stood for “shape files”, but it stands for “simple features”, which actually makes a lot more sense.

Making the tables

I also used another package whose name I thought meant something else in order to make the accompanying summary tables. The tables are easily made with basic dplyr functions, but making them pretty happens with the package {gt}, which stands for “grammar of tables” and not “great tables”, though in this case I think my acronym is superior.

For the simple table in the Quest print, I summarized the number of amenities and their distances and formatted that into a print ready table with {gt}.

# The table covers everything in fun_stuff (the full `fun` list).
fun_sum <- fun_stuff |>
  st_drop_geometry() |>
  group_by(amenity) |>
  summarize(total = n(), avg_dist_miles = mean(dist_miles), sd_dist_miles = sd(dist_miles)) |>
  arrange(desc(total))

# Formatted table, distance as mean ± sd. sd is NA for types with only one place (shows the mean alone).
fun_table <- fun_sum |>
  mutate(amenity = str_to_sentence(str_replace_all(amenity, "_", " "))) |>
  gt() |>
  cols_label(
    amenity = "Type",
    total = "Count",
    avg_dist_miles = "Average Distance (miles)"
  ) |>
  fmt_number(columns = c(avg_dist_miles, sd_dist_miles), decimals = 2) |>
  cols_merge_uncert(col_val = avg_dist_miles, col_uncert = sd_dist_miles) |>  # mean ± sd
  sub_missing(missing_text = "—") |>
  tab_header(
    title = "Fun 'Amenities' Near Reed",
    subtitle = "Distance from the edge of campus"
  ) |>
  tab_source_note("Source: OpenStreetMap contributors") |>
  opt_row_striping() |>
  cols_width(
    amenity ~ px(170),
    total ~ px(110),
    avg_dist_miles ~ px(230)
  ) |>
  cols_align(align = "center", columns = c(total, avg_dist_miles)) |>
  tab_options(data_row.padding = px(10), column_labels.padding = px(10))

For the interactive table, it’s linked to the map using the package {crosstalk}. All you need to do is wrap your dataframe with SharedData$new(), and then any widget built from that shared object responds to the others. So when you pick “ice cream” in the search box, both the map and the table filter at the same time. The table is made with {DT}, which makes sortable, searchable tables. The code here again is a bit easier than the print version, largely because I just accepted the table defaults rather than specifying size and spacing.

# show the name of the thing, the type, distance from Reed, and location 
# (just lat and long, address is likely possible, but I got lazy)
table <- datatable(
  shared,
  colnames = c("Name" = "name", "Type" = "type", "Location (lat, lon)" = "location", "Distance from Reed (mi)" = "dist_miles"),
  rownames = FALSE,
  options = list(
    columnDefs = list(list(targets = c(4, 5), visible = FALSE)),
    order = list(list(3, "asc"))
  )
)

Displaying the interactive map

The last step was saving everything as one HTML file. The {htmlwidgets} package and {htmltools} helped here. Normally, you can just save widgets very easily with the function saveWidget(), but because the graphic is technically several widgets stacked together (the search box, the map, and the table), I had to do another little fiddlesome thing. I wrote a tiny throwaway R Markdown file that renders with “self_contained: true“, so everything ends up in one little file. Then I put that file on GitHub Pages in order to have a spot for it to be hosted so I could drop it into this blog.

Unfortunately, WordPress didn’t want to show it. The plan was to have it embedded on this page with an iframe, which is basically a window into another web page. I did some tinkering to make it display at the correct height so that it could talk to the WordPress page and resize itself accordingly. That works on a normal webpage, but WordPress blogs don’t allow iframes from Github! Only from elite sites like TikTok and Instagram [eyeroll emoji]. But here’s the code that would work if WordPress was slightly better.

# --- Layout and save -------------------------------------------------------
# crosstalk package enables bootstrap styling 
# needed widget to embed into wordpress as an iframe 

# this is supposed to work so the embed can resize itself instead of guessing a
# fixed height and getting an internal scrollbar; to do this the embedding page 
# has to "listen" for this - see amenities/README.md.
resize_reporter <- tags$script(HTML("
  (function() {
    function sendHeight() {
      window.parent.postMessage(
        { type: 'amenities-map-resize', height: document.documentElement.scrollHeight },
        '*'
      );
    }
    window.addEventListener('load', sendHeight);
    new ResizeObserver(sendHeight).observe(document.body);
  })();
"))

widget <- tagList(
  bscols(widths = 12, filter_select("type", "Click in box to display specific amenities", shared, ~type, multiple = TRUE)),
  bscols(widths = 12, legend_box),
  bscols(widths = 12, map),
  bscols(widths = 12, table),
  resize_reporter
)

# this is technically a couple widgets because of the tagList, so can't just use
# saveWidget(), so render through a throwaway .Rmd instead -
# rmarkdown shells out to pandoc, so it needs to be installed
Sys.setenv(RSTUDIO_PANDOC = path.expand("~/.local/bin"))

# this just compiles everything as html with a title and whatnot
rmd_path <- tempfile(fileext = ".Rmd")
writeLines(c(
  "---",
  "title: Ice Cream & Grave Yards - 'Amenities' Near Reed",
  "output:",
  "  html_document:",
  "    self_contained: true",
  "---",
  "",
  "```{r echo=FALSE, warning=FALSE, message=FALSE}",
  "widget",
  "```"
), rmd_path)

# save it to a file
rmarkdown::render(
  rmd_path,
  output_file = "amenities_interactive.html",
  output_dir = "amenities",
  envir = environment(),
  quiet = TRUE
)

Want the full code?

We have a Github repo for this project where you can access everything use to create these maps. It contains two scripts: one for building the image printed in the Quest and one for the interactive map below. From those a map and a table image are created for the Quest and an html file is created for the interactive map.

Open Street Maps has a ton of useful data, so I really encourage people to check it out for a wide variety of projects. And if you ever need help, visit the DataLab!

Looks better in color, huh?

It’s actually surprisingly difficult to convey all the information you’d like to in just black and white and in a small format. I tried grayscale on the map but it didn’t read well, so I just left contour lines, but the colored image shows what’s going on much better.

Chavez Campesinos Blvd

The biggest surprise in making this map was that Chavez has been officially renamed in Oregon’s databases as Campesinos Blvd. They are going to begin changing signs soon. The decision was made official on September 9th by a unanimous city council vote.

How this map was made

To create the contours on this map, I used a file from Oregon Metro’s Regional Land Information System (RLIS). They host hundreds of publicly available datasets that cover a wide range of topics including geographical maps, but files about Oregon demographics, climate and environmental data, civil buildings, transportation, and more. It’s a great resource for many projects.

I used the contour map with five feet intervals to get elevations for the Portland area. This came in the form of shape files (.shp), which initially just looks like a series of lines, not shaded areas between lines.

I wanted an area around Reed, not just Reed itself, so I added a rectangular polygon around the area I wanted. I set it to be white and moved it behind the contour lines.

Because I wanted to include Reed College on my map, I also downloaded another shape file that included names and boundaries for different tax lots in Portland. I opened that as a layer in QGIS and filtered the names to just find the boundary of Reed College and remove all other tax lots. Then I overlaid this Reed layer onto the map.

Then I converted the contour lines into polygons so that they could be shaded in by elevation height. There’s a tool in QGIS called “Contour Polygons” that will do this for you. Then I changed the symbology of that layer so that every five feet difference in elevation appeared as a different color along a gradient. I chose five feet because that is the resolution of the original contour file, but I could have chosen a higher number for a less fine-grained graphic.

Lastly, I wanted to add some roads so that it would be easier to orient oneself. QGIS has a Open Streets Map plugin that allows you to add all kinds of things like bike lanes, horse paths, stop lights, and more. From that layer, I selected the streets that I wanted to add to my map. Then I fiddled with the coloring a bit to make things more legible.

To determine that the height of the Woodstock hill is 240 feet, I used the “identify features” tool to see the information for the contour line that overlapped 39th street.

I did the same for an area on the quad. I identified the quad by adding on a layer of Reed buildings, but I thought the map looked too busy with those included, so I removed them for the final graphic. The point of my graphic was to look at the steepness of the Woodstock hill, so the buildings distracted rather than added to the map.