We’ll first load in variables, which I collected and manipulated in previous scripts.

flow_lines_arc <- st_sf(readRDS("../../results/flow_lines_arc.RDS"))
balt_hex <- st_sf(readRDS("../../results/balt_hex.RDS"))
otm_lehd <- st_read("../../results/otm_5b749ca22bf64d9abe899ef50a619131/points_2019.shp") %>% st_transform(3857) %>%
  select(id, c000, geometry)
scooters_raw <- st_read("../../results/scooters_raw_050122_050722.gpkg") %>% filter(is_disabled == 0)

Now, we’ll use tidycensus and tigris to grab census data on race, ethnicity, income, and housing units at the block group level. We’ll also get geometries for Baltimore City and County for visualization.

baltimore_bg_income <- get_acs(geography = "block group", 
                           variables = c("pop" = "B03002_001", # Total
                                         "pop_nhwhite" = "B03002_003", # NH White
                                         "pop_nhblack" = "B03002_004", # NH Black
                                         "pop_nhamind" = "B03002_005", # NH Am Ind
                                         "pop_nhasian" = "B03002_006", # NH Asian
                                         "pop_nhhwnpi" = "B03002_007", # NH Hawaiin/PI
                                         "pop_nhother" = "B03002_008", # One Other
                                         "pop_nhtwomr" = "B03002_009", # Two+
                                         "pop_hispltx" = "B03002_012", # Hispanic/Latinx
                                         "hu_total"  = "B25001_001", # Housing Units
                                         "hu_totocc" = "B25003_001", # Housing Units - Occ
                                         "hu_totown" = "B25003_002", # Housing Units - Owner Occ,
                                         "hu_totrnt" = "B25003_003", # Housing Units - Renter Occ,
                                         "mhhi" = "B19013_001"), #median household income
                           
                           year = 2019,
                           survey = "acs5",
                           state = c(24), 
                           county = c(510, 005), 
                           geometry = TRUE, 
                           output = "wide",
                           progress_bar = FALSE) %>% st_transform(3857)
balt_city_bound <- counties("MD", cb = TRUE, progress_bar = FALSE) %>% filter(COUNTYFP == 510) %>% st_transform(3857)

We now conflate race and ethnicity, though problematic, to make calculations more simple.

# Computes the NH Asian Population
baltimore_bg_income$pop_nhasianXE <- baltimore_bg_income$pop_nhasianE + baltimore_bg_income$pop_nhhwnpiE

# Computes the NH "Other" Population
baltimore_bg_income$pop_nhotherXE <- baltimore_bg_income$pop_nhamindE + baltimore_bg_income$pop_nhotherE + baltimore_bg_income$pop_nhtwomrE

We define a function, st_erase, to remove water area (again for visualization purposes).

st_erase <- function(x, y) {
  st_difference(x, st_make_valid(st_union(y)))}
bmore_water <- area_water("MD", c(510,005), class = "sf") %>% st_transform(3857) %>% filter(AWATER > 20000)
bmore_water <- st_make_valid(st_buffer(bmore_water, 0)) # Fix topology

We use our st_erase function to erase water from Baltimore City’s boundary, then begin our areal interpolation (hexbinning) operations. Since our LEHD data is represented as points, we can perform a spatial join and then sum the jobs within each hex. In this code block, we also define a subset of our hex grid, which we’ll use to visualize only the city in our maps.

baltimore_bg_income <- baltimore_bg_income %>% st_transform(3857)
balt_city_bound <- st_erase(balt_city_bound, bmore_water)
balt_hex <- balt_hex %>% st_transform(3857)
balt_bg_emp <- st_join(otm_lehd, balt_hex)
balt_bg_emp <- balt_bg_emp %>% group_by(hex_id) %>% summarise(sum_c000 = sum(c000))
balt_hex.intersects <- st_intersects(st_union(balt_city_bound), balt_hex)
balt_hex.subset <- balt_hex[balt_hex.intersects[[1]],]
bmore_scooters_in_hex <- st_join(balt_hex.subset, scooters_raw)
bmore_scooter_hex_count <- count(as_tibble(bmore_scooters_in_hex), hex_id)
balt_hex.subset <- left_join(balt_hex.subset, bmore_scooter_hex_count)

Since our population data is represented as areal objects (block group polygons), we can’t perform a simple spatial join: we must use the areal package to interpolate the population per hex.

People-points

I’ve coined a new unit: people-points. Represented by one job as represented in LEHD data, added with the total population in a given area, people-points can be used to find centers of transport and human activity as people commute to and from work and home. I’ll be symbolizing my maps based on the number of scooters per person-point (in this case, per 1,000 people-points for better scaling) in order to pick out locations where the number of scooters is not proportional to the number of jobs and residents in an area.

ar_validate(source = baltimore_bg_income, target = balt_hex.subset, varList = "popE", method = "aw", verbose = TRUE)

bmore_grid_income <- aw_interpolate(balt_hex.subset, tid = hex_id, source = baltimore_bg_income, sid = "GEOID", weight = "sum", output = "sf", extensive = c("popE","hu_totalE", "pop_nhwhiteE"))
bmore_grid_emp_income <- left_join(bmore_grid_income, st_drop_geometry(balt_bg_emp))
bmore_grid_emp_income[is.na(bmore_grid_emp_income)] <- 0
bmore_grid_emp_income <- bmore_grid_emp_income %>% 
                         rename(scooter_count = n)%>% 
                         mutate(people_pt = popE + sum_c000, 
                                scooters_per_pt = (1000*scooter_count/126)/people_pt, #divide by number of 15-min intervals in time collected
                                pct_bipoc = (popE-pop_nhwhiteE)/popE) %>% 
                                filter(people_pt > 100)
bmore_grid_emp_income[is.na(bmore_grid_emp_income)] <- 0
bmore_grid_emp_income <- bmore_grid_emp_income %>% st_transform(3857)

We’ll set our bounding box so only Baltimore City is visible in the maps we create…

bbox_new <- st_bbox(bmore_grid_emp_income) # current bounding box

xrange <- bbox_new$xmax - bbox_new$xmin # range of x values
yrange <- bbox_new$ymax - bbox_new$ymin # range of y values

bbox_new <- bbox_new %>%  # take the bounding box ...
  st_as_sfc() # ... and make it a sf polygon

And then we’ll create our maps!

Maps

I used ggplot2 to create these maps, though I generally prefer tmap, so I could use patchwork for easy arrangement.

layout1 <- ggplot.scooters_pt|ggplot.pct_bipoc
layout1

layout2 <- ggplot.scooters_pt | ggplot.pts
layout2

We can see that, during a week of morning commuting (Monday - Monday), the areas with the most scooter trips per person-point are in Locust Point and the Inner Harbor. More specifically, the two yellowest polygons, which represent the hexes with the highest proportion of trips per person-point, contain the Under Armour main campus. Notably, areas with high percentages of BIPOC individuals are almost all dark purple, showing few numbers of scooter trips per person-point.

What does this mean?

According to Spin, they are “committed to being the best possible partner for cities while building the safest, most equitable, and most sustainable mobility solution for the communities [they] serve” (emphasis added). I’d say that, whether intentional or not, the stark contrast of micromobility in minority communities compared to whiter communities demonstrates a veritable lack of equitable access, at least in Baltimore. Whether the vehicles are purposely dropped in whiter areas after charging, if they’re used to commute from these areas to Downtown but not back, or any other reason, calling the distribution of scooter trips and availability “equitable” is downright laughable as things stand at the moment.

Further research:

An important reason that cities (including Baltimore) adopt these platforms is to supplement existing transit networks by providing a method for a individual to quickly get to a transit stop before switching to the bus, train, or other method of transit. One interesting area of research in this regard could be to examine how many scooter trips have start and end points that mirror existing transit routes. Could people be using scooters instead of public transit? If so, why? And, as always, more research could be carried out with more data. Since the script collecting our data runs every 15 minutes until the server runs out of storage or the inevitable collapse of capitalism (whichever comes first), we could theoretically run this same analysis for time periods of months or even years, given powerful enough hardware.