Lab 8 Deliverable
library(tidyverse)
library(sf)
library(sp)
library(tmap)
library(osmdata)
library(ggmap)
library(readxl)
library(classInt)
library(cowplot)
library(maps)
library(ggspatial)
library(basemaps)
library(mapview)Data for the Maps
tree_shapes <- st_read("data/Parks_Tree_Inventory-shp/")## Reading layer `fd3ad191-4bf0-4577-b6ec-724baeb619792020413-1-a3eocm.isio' from data source `C:\Users\boyernat\OneDrive - Oregon Health & Science University\Documents\2025-04 Spring Term\BMI 525 Data Visualization\Labs\08_lab\data\Parks_Tree_Inventory-shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 25534 features and 39 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -13667450 ymin: 5692149 xmax: -13635010 ymax: 5724486
## Projected CRS: WGS 84 / Pseudo-Mercatorpdx_boundaries <- st_read("data/Neighborhoods__Regions_-shp/")## Reading layer `Neighborhoods__Regions_' from data source
## `C:\Users\boyernat\OneDrive - Oregon Health & Science University\Documents\2025-04 Spring Term\BMI 525 Data Visualization\Labs\08_lab\data\Neighborhoods__Regions_-shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 98 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -13677560 ymin: 5689875 xmax: -13632920 ymax: 5724919
## Projected CRS: WGS 84 / Pseudo-Mercatorriver_boundaries <- st_read("data/Willamette_Columbia_River_Ordinary_High_Water-shp/")## Reading layer `Willamette_Columbia_River_Ordinary_High_Water' from data source
## `C:\Users\boyernat\OneDrive - Oregon Health & Science University\Documents\2025-04 Spring Term\BMI 525 Data Visualization\Labs\08_lab\data\Willamette_Columbia_River_Ordinary_High_Water-shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -13676470 ymin: 5686111 xmax: -13632810 ymax: 5727803
## Projected CRS: WGS 84 / Pseudo-MercatorChallenge 1
Map centered on your neighborhood!
pdx_boundaries %>%
ggplot() +
geom_sf() +
coord_sf(
xlim=c(-13658000,-13650000),
ylim=c(5699500, 5707500)) +
geom_sf_label(aes(label = MAPLABEL), size = 3) +
geom_sf_label(data = pdx_boundaries %>% filter(MAPLABEL == "Buckman Community Association"),
aes(label = MAPLABEL,
fill = "pink"),
size = 3) +
theme(legend.position = "none")
Challenge 2
Changing river and neighborhood color!
First, simplify river points for quicker plotting.
river_boundaries <- river_boundaries %>% st_simplify(dTolerance = 10)pdx_boundaries %>%
ggplot() +
geom_sf(fill = "palegreen1",
col = "black") +
geom_sf(data = river_boundaries,
fill = "lightblue1",
col = "midnightblue",
size = 0.0) +
theme_classic()
Challenge 3
Experiment with projections!
# world map
world <- st_as_sf(map('world', plot = FALSE, fill = TRUE))New projection: “August Epicycloidal” Attempting to center on the US – can see that this projection creates some major distortion for areas on opposite side of the globe from the US.
world2 <- st_transform(
world,
"+proj=august +ellps=WGS84 +over +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0"
)
ggplot() +
geom_sf(data = world2)## Warning in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy, :
## GDAL Error 1: PROJ: pipeline: Pipeline: Inverse operation for august is not
## available
Now, the coordinates are flipped (opposite sign) in an attempt to focus of the other side of the world. A similar problem persists.
world3 <- st_transform(
world,
"+proj=august +ellps=WGS84 +over +lat_1=-29.5 +lat_2=-45.5 +lat_0=-37.5 +lon_0=96 +x_0=0 +y_0=0"
)
ggplot() +
geom_sf(data = world3)## Warning in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy, :
## GDAL Error 1: PROJ: pipeline: Pipeline: Inverse operation for august is not
## available
Challenge 4
Plot both grocery stores and farmers markets!
Data:
farmers_market <- st_read("data/Farmers_Markets-shp/")## Reading layer `Farmers_Markets' from data source
## `C:\Users\boyernat\OneDrive - Oregon Health & Science University\Documents\2025-04 Spring Term\BMI 525 Data Visualization\Labs\08_lab\data\Farmers_Markets-shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 20 features and 11 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -13665000 ymin: 5696460 xmax: -13644270 ymax: 5714893
## Projected CRS: WGS 84 / Pseudo-Mercatorgrocery_stores <- st_read("data/Grocery_Stores-shp/")## Reading layer `Grocery_Stores' from data source
## `C:\Users\boyernat\OneDrive - Oregon Health & Science University\Documents\2025-04 Spring Term\BMI 525 Data Visualization\Labs\08_lab\data\Grocery_Stores-shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 123 features and 9 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -13665520 ymin: 5686206 xmax: -13634020 ymax: 5718088
## Projected CRS: WGS 84 / Pseudo-MercatorPlot:
pdx_boundaries %>%
ggplot() +
geom_sf(fill = "white", col = 'black') +
geom_sf(data = river_boundaries,
fill = "lightblue1", size = 0.0, color = 'midnightblue') +
geom_sf(data = farmers_market,
mapping = aes(color = Day),
shape = 8) +
labs(color = "Farmer's Markets Days Open") +
scale_color_brewer(palette = "Dark2") +
ggnewscale::new_scale_color() +
geom_sf(data = grocery_stores,
mapping = aes(color = TYPE, shape = TYPE)) +
labs(shape = "Grocery Store Type",
color = "Grocery Store Type") +
scale_color_manual(values = c('Large Chain Grocery' = 'darkblue', 'Independent or Ethnic Grocery' = 'red3')) +
theme_classic()
Challenge 5
Plot your neighborhood’s park
Getting park coordinates with drawing box help.
laurelhurst_shape <- basemaps::draw_ext()## Loading required namespace: leaflet.extras## Loading required package: shiny##
## Listening on http://127.0.0.1:4466# check basemap's coordinate system
st_crs(laurelhurst_shape)## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## MEMBER["World Geodetic System 1984 (G2296)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]# check pdx's coordinate system
st_crs(pdx_boundaries)## Coordinate Reference System:
## User input: WGS 84 / Pseudo-Mercator
## wkt:
## PROJCRS["WGS 84 / Pseudo-Mercator",
## BASEGEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## MEMBER["World Geodetic System 1984 (G2296)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]],
## CONVERSION["Popular Visualisation Pseudo-Mercator",
## METHOD["Popular Visualisation Pseudo Mercator",
## ID["EPSG",1024]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["False easting",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["easting (X)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["northing (Y)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Web mapping and visualisation."],
## AREA["World between 85.06°S and 85.06°N."],
## BBOX[-85.06,-180,85.06,180]],
## ID["EPSG",3857]]# transform laurelhurst coords to pdx's coord
laurelhurst_tranformed <- st_transform(laurelhurst_shape, st_crs(pdx_boundaries))
# get bounding box from laurelhurst
laurelhurst_coords <- st_bbox(laurelhurst_tranformed)
# bounding box with default stadia map coords
laurelhurst_coords2 <- st_bbox(st_transform(laurelhurst_shape, crs = 4326))Plotting Laurelhurst Par and its trees:
pdx_boundaries %>% ggplot() +
geom_sf() +
geom_sf(data = tree_shapes,
mapping = aes(size = CrownWidth,
color = as.numeric(as.factor(Species))),
shape = 8) +
scale_color_gradientn(colors = c("#d9f0a3", "#78c679", "#238443", "#004529")) +
guides(color = FALSE,
size = FALSE) +
scale_size(range = c(0.05, 2)) +
coord_sf(xlim = laurelhurst_coords[c("xmin", "xmax")],
ylim = laurelhurst_coords[c("ymin", "ymax")]) +
theme_classic() +
labs(title = "Laurelhurst Park Tree Inventory",
subtitle = "Trees colored by species, size indicative of crown width") 
Challenge 6
Use OpenStreetMap tiles
# stadia map on Laurelhurst park
laurelhurst_basemap <- get_stadiamap(c(left=-122.63245, bottom=45.51824, right=-122.62059, top=45.52357),
zoom = 16,
maptype="stamen_terrain")## ℹ © Stadia Maps © Stamen Design © OpenMapTiles © OpenStreetMap contributors.# Tree map coords tranformed for stadia
trees_transformed <- st_transform(tree_shapes, st_crs(4326))
ggmap(laurelhurst_basemap) +
geom_sf(data = trees_transformed,
inherit.aes = FALSE,
mapping = aes(size = CrownWidth,
color = as.numeric(as.factor(Species))),
shape = 8) +
scale_color_gradientn(colors = c("#d9f0a3", "#78c679", "#238443", "#004529")) +
guides(color = FALSE,
size = FALSE) +
scale_size(range = c(0.05, 2)) +
theme_classic() +
labs(title = "Laurelhurst Park Tree Inventory",
subtitle = "Trees colored by species, size indicative of crown width",
y = "Latitude",
x = "Longitude") ## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.## Warning: Removed 263 rows containing missing values or values outside the scale range
## (`geom_sf()`).
Challenge 7
A map using a different census column (or multiple columns!)
# read in census data, selecting up until desired column and rows
pdx_population <- read_excel("data/Census_2010_Data_Cleanedup.xlsx",
sheet="Census_2010_Neighborhoods",
range="A7:CM101",
col_names = FALSE)
# columns to keep
c1 <- 1
c2 <- ncol(pdx_population) - 3
c3 <- ncol(pdx_population)
keep_cols <- c(c1, c2, c3)
# subset columns of interest (first = neighborhoods, last = renters)
pdx_population <- pdx_population[, keep_cols] %>%
rename("Neighborhood" = "...1", "rentals" = "...91", "total_homes" = "...88") %>%
mutate(Neighborhood=as.factor(Neighborhood)) %>%
mutate(rental_density = rentals / total_homes)
# we need to do a tiny bit of cleanup; the census spreadsheet has some slightly different names
pdx_population <- pdx_population %>% mutate(Neighborhood=recode(Neighborhood,
"ARGAY" = "ARGAY TERRACE",
"BROOKLYN" = "BROOKLYN ACTION CORPS",
"BUCKMAN" = "BUCKMAN COMMUNITY ASSOCIATION",
"CENTENNIAL" = "CENTENNIAL COMMUNITY ASSOCIATION",
"CULLY" = "CULLY ASSOCIATION OF NEIGHBORS",
"CENTENNIAL" = "CENTENNIAL COMMUNITY ASSOCIATION",
"DOWNTOWN" = "PORTLAND DOWNTOWN",
"GOOSE HOLLOW" = "GOOSE HOLLOW FOOTHILLS LEAGUE",
"HAYDEN ISLAND" = "HAYDEN ISLAND NEIGHBORHOOD NETWORK",
"HOSFORD-ABERNETHY" = "HOSFORD-ABERNETHY NEIGHBORHOOD DISTRICT ASSN.",
"IRVINGTON" = "IRVINGTON COMMUNITY ASSOCIATION",
"LLOYD DISTRICT" = "LLOYD DISTRICT COMMUNITY ASSOCIATION",
"NORTHWEST DISTRICT" = "NORTHWEST DISTRICT ASSOCIATION",
"OLD TOWN-CHINATOWN" = "OLD TOWN COMMUNITY ASSOCIATION",
"PARKROSE HEIGHTS" = "PARKROSE HEIGHTS ASSOCIATION OF NEIGHBORS",
"PEARL" = "PEARL DISTRICT",
"SABIN" = "SABIN COMMUNITY ASSOCIATION",
"SELLWOOD-MORELAND" = "SELLWOOD-MORELAND IMPROVEMENT LEAGUE",
"SOUTHWEST HILLS" = "SOUTHWEST HILLS RESIDENTIAL LEAGUE",
"SUMNER" = "SUMNER ASSOCIATION OF NEIGHBORS",
"SUNDERLAND" = "SUNDERLAND ASSOCIATION OF NEIGHBORS",
"WILKES" = "WILKES COMMUNITY GROUP"
))
# and now left join with our boundary data set- note that there will be some regions with no population data, so we want a left-join
boundaries_with_pop <- left_join(pdx_boundaries, pdx_population, by=c("NAME"="Neighborhood"))
boundaries_with_pop %>% glimpse## Rows: 98
## Columns: 14
## $ OBJECTID <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, …
## $ NAME <chr> "CATHEDRAL PARK", "UNIVERSITY PARK", "PIEDMONT", "WOODL…
## $ COMMPLAN <chr> NA, NA, "ALBINA", "ALBINA", NA, "ALBINA", "ALBINA", "AL…
## $ SHARED <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ COALIT <chr> "NPNS", "NPNS", "NPNS", "NECN", "CNN", "NPNS", "NPNS", …
## $ HORZ_VERT <chr> "HORZ", "HORZ", "VERT", "HORZ", "HORZ", "HORZ", "HORZ",…
## $ MAPLABEL <chr> "Cathedral Park", "University Park", "Piedmont", "Woodl…
## $ ID <int> 31, 88, 70, 93, 23, 2, 66, 19, 67, 84, 44, 48, 89, 63, …
## $ Shape_Leng <dbl> 11434.255, 11950.860, 10849.327, 8078.361, 18179.392, 9…
## $ Shape_Area <dbl> 5424297.8, 6981456.8, 6079530.0, 3870553.6, 16580624.7,…
## $ total_homes <dbl> 1703, 1834, 2983, 2084, 5230, 2774, 2756, 4001, 2314, 8…
## $ rentals <dbl> 699, 617, 1030, 656, 2148, 825, 919, 1125, 1160, 222, 1…
## $ rental_density <dbl> 0.41045214, 0.33642312, 0.34528998, 0.31477927, 0.41070…
## $ geometry <MULTIPOLYGON [m]> MULTIPOLYGON (((-13664216 5..., MULTIPOLYG…boundaries_with_pop %>%
ggplot() +
geom_sf(aes(fill = rental_density)) +
geom_sf(data = river_boundaries, fill = 'lightblue1', size = 0) +
theme_classic() +
labs(title = "Portland, OR Neighborhood Rentals",
subtitle = "Proportion of homes that are occupied by renters, by neighborhood",
y = "Latitude",
x = "Longitude",
fill = "Proportion of Rentals")
Challenge 8
Exploring color scales
Distribution of rental proportions.
ggplot(data = boundaries_with_pop,
mapping = aes(x = reorder(NAME, -rental_density),
y= rental_density)) +
geom_col(col = "white") +
theme_classic() +
theme(axis.text.x = element_text(angle = 90,
vjust = 0.5,
hjust = 1,
size = 3)) +
labs(y = "Proportion of Rentals",
x = "Neighborhoods")
ggplot(data = boundaries_with_pop,
mapping = aes(x = rental_density)) +
geom_histogram() +
labs(x = "Proprotion of Rentals")
Different approaches to discretizing rental density.
# set up
min.rental <- min(boundaries_with_pop$rental_density, na.rm = TRUE)
max.rental <- max(boundaries_with_pop$rental_density, na.rm = TRUE)
range.rental <- max.rental - min.rental
stdev.rental <- sd(boundaries_with_pop$rental_density, na.rm = TRUE)
# Equal interval bins (bins = 7)
equal.interval <- seq(min.rental, max.rental, by = range.rental/6)
# Quantile interval bins
quantile.interval <- quantile(boundaries_with_pop$rental_density, probs=seq(0,1, by = 1/7), na.rm = TRUE)
# Standard deviation intervals
std.interval <- c(seq(min.rental, max.rental, by = stdev.rental), max.rental)
# "Jenk's" Intervals
jenks.intervals <- classIntervals(boundaries_with_pop$rental_density, n=7, style = 'jenks')$brks## Warning in classIntervals(boundaries_with_pop$rental_density, n = 7, style =
## "jenks"): var has missing values, omitted in finding classes# Assign breaks to variable
# Equal interval
boundaries_with_pop$rental.equal <- cut(boundaries_with_pop$rental_density, breaks = equal.interval, include.lowest = TRUE)
# Quantile intervals
boundaries_with_pop$rental.quantile <- cut(boundaries_with_pop$rental_density, breaks = quantile.interval, include.lowest = TRUE)
# St Dev Intervals
boundaries_with_pop$rental.sd <- cut(boundaries_with_pop$rental_density, breaks = std.interval, include.lowest = TRUE)
# Jenk's intervals
boundaries_with_pop$rental.jenks <- cut(boundaries_with_pop$rental_density, breaks = jenks.intervals, include.lowest = TRUE)Plotting these different intervals
barcharts <- function(break_col) {
boundaries_with_pop %>%
filter(!is.na(rental_density)) %>%
ggplot(mapping = aes(x = fct_reorder(MAPLABEL, -rental_density),
y = rental_density)) +
geom_col(aes(fill = .data[[break_col]])) +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank()) +
scale_y_continuous(name = "Proportion of Rentals") +
scale_x_discrete(name = NULL) +
scale_fill_discrete(guide = FALSE) +
ggtitle(break_col)
}
plot_grid(
barcharts("rental.equal"),
barcharts("rental.quantile"),
barcharts("rental.sd"),
barcharts("rental.jenks"),
nrow = 2,
ncol = 2
)## Warning: The `guide` argument in `scale_*()` cannot be `FALSE`. This was deprecated in
## ggplot2 3.3.4.
## ℹ Please use "none" instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Applying these intervals to the map.
map_intervals <- function(break_col) {
boundaries_with_pop %>%
ggplot() +
geom_sf(aes(fill = .data[[break_col]])) +
geom_sf(data = river_boundaries, fill = 'lightblue1', size = 0) +
theme_classic() +
labs(title = break_col,
subtitle = "Proportion of homes that are occupied by renters, by neighborhood",
y = "Latitude",
x = "Longitude",
fill = "Proportion of Rentals") +
scale_fill_brewer(type = 'seq',
palette = 'RdPu')
}
map_intervals("rental.equal")
map_intervals("rental.quantile")
map_intervals("rental.sd")
map_intervals("rental.jenks")