The Cackling Goose and the Canada Goose: can you tell them apart if they aren’t right next to each other? People at Oaks Bottom probably can’t.

This week we looked at the geese present on and around Reed’s campus. Anecdotally, I think there are more geese on campus than there were when I was here in the mid to late 2000s. If we want to look at confirming or denying that assertion there are a number of routes that we could take. I’m going to go through one that I think might be people’s first thought for addressing this question and point out some potential pitfalls in what might seem like an otherwise straightforward analysis.

Some of the problems discussed here will be specific to the topic, but many can be broadly applied to other types of research. When you’re working on a project, it’s always good to check your assumptions to make sure you thoroughly understand what’s in your data, including how it is collected, potential confounding variables, and issues with masking or silent imputation (for example, how null data is handled).

Public Data & Citizen Science

Public, online repositories of data are extremely useful for all kinds of research. OpenStreetMaps has an extensive database for maps and location labeling of all kinds (as described a little in the previous blog entry). PurpleAir and OpenAQ have data on air quality that is fine-scaled and up to the minute because of user-submitted measurements. Globe at Night is a citizen science initiative to measure light pollution around the world and relies totally on individual contributions. These types of organizations provide a wealth of useful knowledge that is available for anyone to use, but you have to be careful how you interpret it.

This project relied on eBird data. eBird is a site for birders to submit lists of birds that they saw or heard, how many there were, where they were, and a number of other things. Birders are often obsessive nutjobs¹, so that makes eBird one of, if not the, the largest repository of biodiversity data in existence. Scientists use it to track migration, population trends, and the effects of climate change and habitat loss, and anyone can explore the maps and species data online. You can look at birding hotspots worldwide or right in your neighborhood, like the one at Reed Lake (aka, the canyon).

eBird is a bit cumbersome to use, so I accessed its data via the Global Biodiversity Information Facility, which pulls data from a number of open access sources (it’s very cool, check it out). To do this I used an R package called {rgbif} to pull data from three eBird hotspots near Reed (Reed Lake, Oaks Bottom, and the Rhodie Gardens) for both Cackling Geese (Branta hutchinsii) and Canada Geese (Branta canadensis). That code is in the ‘pull_geese.R’ script on the GitHub, but the gist of it relies on having a free user key for eBird and this function:

occ_search(
    taxonKey   = key,
    datasetKey = ebird_key,
    locality   = loc_id,
    year       = year_range,
    limit      = 20000
  )

This pulls data since 2010 for the three locations (as identified by unique eBird IDs), and that data comes as individual people’s checklists, so a single row, looks something like this and includes the specific date, whose list it is, and how many birds they saw.

By summing up that data by year, you can get the following figure, which is similar to what is in The Quest, but I added the Rhodie Gardens just for fun.

On the face of it, this is a simple and straightforward way to look at the number of geese over time. But there are a number of issues, some obvious and some not so obvious, that make this graph not as useful as one might initially think.

Problems

The easiest data to use for analysis is highly structured data. This is data that has a design implemented before collection, has a fixed scheme for how it is recorded, and has neat rows and columns. eBird data may appear this way at first but it is actually semi-structured data. The form for entering the bird count data might be the same for everyone, but participants get to choose where, when, and for how long they want to record. They can log every bird they see or just the ones they’re interested in. And their observations are based solely on their own expertise.

Structured vs Semi-Structured Data

An example of similar data that is structured is the USGS North American Breeding Bird Survey. Data is still collected by volunteers, but those volunteers have to officially register with their local government and follow a strict protocol. The routes are random, the survey takes place at the same time each year, participants count for a set length of time, and the skill levels are roughly equally advanced. Because of this every count represents a standardized effort, so year-to-year differences reflect real changes in bird populations rather than differences in how the data was collected.

Let’s dig into the issues.

Sampling Bias

Sampling bias occurs when a sample or survey is not representative of the whole population. As you might expect, people are going to bird in nice, easy to reach places rather than rural, inhospitable, or just plain ugly locations. So eBird checklists cluster near roads, cities, parks, and hotspots, and around wealthy regions where birds are abundant. Just looking at a map of where records are, you can see obvious oversampling in some areas and no data for others. Here’s a map of an area around Portland.

An absence of records in a location obviously doesn’t mean there are no birds there. We expect to see higher numbers of birds recorded in heavily trafficked areas. So, for our goose data, do we really think that geese flock to the Rhodie Gardens far more than to Reed? Perhaps, but think about how likely the average person on Reed campus is to be recording bird totals versus the average visitor to the Rhodie Gardens.

Continue reading →

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.