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(ggthemes)
library(leaflet)
library(leaflet.esri)

Challenge 1

Map centered on your neighborhood!

pdx_boundaries <- st_read("data/Neighborhoods__Regions_-shp")

## Reading layer `Neighborhoods__Regions_' from data source `/cloud/project/data/Neighborhoods__Regions_-shp' using driver `ESRI Shapefile'
## Simple feature collection with 98 features and 10 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -13677560 ymin: 5689875 xmax: -13632920 ymax: 5724919
## CRS:            3857

river_boundaries <- st_read("data/Willamette_Columbia_River_Ordinary_High_Water-shp/")

## Reading layer `Willamette_Columbia_River_Ordinary_High_Water' from data source `/cloud/project/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
## bbox:           xmin: -13676470 ymin: 5686111 xmax: -13632810 ymax: 5727803
## CRS:            3857

Choosing Mulnomah Neighborhood

Multnomah <- pdx_boundaries %>% arrange(-Shape_Area) %>% filter(MAPLABEL=='Multnomah')

pdx_boundaries %>% ggplot() + geom_sf() + 
  coord_sf(
    xlim=c( -13665229,-13657000), 
    ylim=c(5690426, 5698669)) +
    geom_sf_label(aes(label=MAPLABEL))

Challenge 2

Changing river and neighborhood color!

I was once told Wilammette means “grean snake” in the language of native americans, who lived here and named the river for it’s color. I don’t know if it is at all credible, but that is why I chose the turquoise color for the river.

pdx_boundaries %>% ggplot() + geom_sf(fill='navajowhite3') + 
    geom_sf(data=river_boundaries, fill="turquoise3", size=0.0)+
  coord_sf(
         xlim=c( -13664656,-13652123), 
         ylim=c(5689075, 5699953))+
  theme_economist()+
  geom_sf_label(data=Multnomah, aes(label=MAPLABEL))

Challenge 3

Experiment with projections!

Looking at the world with my native city. Saint Petersburg, Russia is in the center:

world1 <- sf::st_as_sf(map('world', plot = FALSE, fill = TRUE))
world2 <- sf::st_transform(
  world1,
  "+proj=laea +y_0=0 +lon_0=30.3609 +lat_0=59.9311  +ellps=WGS84 +no_defs"
)
ggplot() + geom_sf(data = world2)

Here is looking at another Portland, Main

world2 <- sf::st_transform(
  world1,
  "+proj=laea +y_0=0 +lon_0=-70.2568 +lat_0=43.659  +ellps=WGS84 +no_defs"
)
ggplot() + geom_sf(data = world2)

Challenge 4

Plot both grocery stores and farmers markets!

trees_shape <- st_read("data/Parks_Tree_Inventory-shp")

## Reading layer `fd3ad191-4bf0-4577-b6ec-724baeb619792020413-1-a3eocm.isio' from data source `/cloud/project/data/Parks_Tree_Inventory-shp' using driver `ESRI Shapefile'
## Simple feature collection with 25534 features and 39 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -13667450 ymin: 5692149 xmax: -13635010 ymax: 5724486
## CRS:            3857

farmers_market <- st_read("data/Farmers_Markets-shp/")

## Reading layer `Farmers_Markets' from data source `/cloud/project/data/Farmers_Markets-shp' using driver `ESRI Shapefile'
## Simple feature collection with 20 features and 11 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -13665000 ymin: 5696460 xmax: -13644270 ymax: 5714893
## CRS:            3857

grocery_stores <- st_read("data/Grocery_Stores-shp/")

## Reading layer `Grocery_Stores' from data source `/cloud/project/data/Grocery_Stores-shp' using driver `ESRI Shapefile'
## Simple feature collection with 123 features and 9 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -13665520 ymin: 5686206 xmax: -13634020 ymax: 5718088
## CRS:            3857

pdx_boundaries %>% ggplot() + geom_sf(fill='palegoldenrod') + 
    geom_sf(data=river_boundaries, fill="turquoise3", size=0.0)+
    geom_sf(data=farmers_market, mapping=aes(color=Day))+
  coord_sf(
         xlim=c( -13664656,-13652123), 
         ylim=c(5689075, 5699953))+
  geom_sf_label(data=Multnomah, aes(label=MAPLABEL))+
  ggtitle("Farmer's Markets near Multnomah by Day")+
  theme_economist()

pdx_boundaries %>% ggplot() +
  geom_sf(fill='palegoldenrod') + 
  geom_sf(data=river_boundaries, fill="turquoise3", size=0.0)+
  geom_sf(data=grocery_stores, mapping=aes(shape=STATUS, color=TYPE ))+
  coord_sf(
         xlim=c( -13664656,-13652123), 
         ylim=c(5689075, 5699953))+
  geom_sf_label(data=Multnomah, aes(label=MAPLABEL))+
  ggtitle("Grocery Stores near Multnomah by Type and Status")+
  theme_economist()

Challenge 5

Plot your neighborhood’s park

I’ve chosen Wilamette Park

nrow(trees_shape)

## [1] 25534

pdx_boundaries %>% ggplot() + geom_sf() + 
  geom_sf(data=river_boundaries, fill="blue", size=0.0) + 
  geom_sf(data=trees_shape)

#Exploring trees

trees_raw <- read_csv("data/Parks_Tree_Inventory.csv")

## Parsed with column specification:
## cols(
##   .default = col_character(),
##   X = col_double(),
##   Y = col_double(),
##   OBJECTID = col_double(),
##   DBH = col_double(),
##   TreeHeight = col_double(),
##   CrownWidthNS = col_double(),
##   CrownWidthEW = col_double(),
##   CrownBaseHeight = col_double(),
##   UserID = col_double(),
##   Structural_Value = col_double(),
##   Carbon_Storage_lb = col_double(),
##   Carbon_Storage_value = col_double(),
##   Carbon_Sequestration_lb = col_double(),
##   Carbon_Sequestration_value = col_double(),
##   Stormwater_ft = col_double(),
##   Stormwater_value = col_double(),
##   Pollution_Removal_value = col_double(),
##   Pollution_Removal_oz = col_double(),
##   Total_Annual_Benefits = col_double()
## )

## See spec(...) for full column specifications.

glimpse(trees_raw)

## Observations: 25,534
## Variables: 41
## $ X                          <dbl> -13658194, -13658213, -13658254, -13658227…
## $ Y                          <dbl> 5712487, 5712485, 5712491, 5712487, 571248…
## $ OBJECTID                   <dbl> 426, 427, 428, 429, 430, 431, 432, 433, 43…
## $ Inventory_Date             <chr> "2017/05/09 00:00:00+00", "2017/05/09 00:0…
## $ Species                    <chr> "PSME", "PSME", "CRLA", "QURU", "PSME", "P…
## $ DBH                        <dbl> 37.4, 32.5, 9.7, 10.3, 33.2, 32.1, 28.4, 2…
## $ Condition                  <chr> "Fair", "Fair", "Fair", "Poor", "Fair", "F…
## $ TreeHeight                 <dbl> 105, 94, 23, 28, 102, 95, 103, 105, 97, 11…
## $ CrownWidthNS               <dbl> 44, 49, 28, 38, 43, 35, 40, 45, 56, 35, 45…
## $ CrownWidthEW               <dbl> 57, 45, 27, 31, 44, 39, 40, 29, 45, 33, 50…
## $ CollectedBy                <chr> "staff", "staff", "staff", "staff", "staff…
## $ Notes                      <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ CrownBaseHeight            <dbl> 4, 4, 3, 5, 4, 12, 6, 13, 5, 17, 6, 6, 4, …
## $ GlobalID                   <chr> "{8413C6D7-038A-4685-AFCB-C3BCE4D90A24}", …
## $ StaffReview                <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ created_user               <chr> "PARKS_ADMIN", "PARKS_ADMIN", "PARKS_ADMIN…
## $ created_date               <chr> "2017/05/23 19:37:43+00", "2017/05/23 19:3…
## $ last_edited_user           <chr> "PKJOSHUA", "PKJOSHUA", "PKJOSHUA", "PKJOS…
## $ last_edited_date           <chr> "2019/12/11 00:42:17+00", "2019/12/11 00:4…
## $ UserID                     <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,…
## $ Structural_Value           <dbl> 10101.06, 7900.55, 1110.81, 525.26, 8224.5…
## $ Carbon_Storage_lb          <dbl> 3223.3, 2460.8, 310.3, 400.1, 2533.4, 2333…
## $ Carbon_Storage_value       <dbl> 209.07, 159.62, 20.13, 25.95, 164.33, 151.…
## $ Carbon_Sequestration_lb    <dbl> 36.8, 31.3, 14.8, 6.0, 32.0, 30.8, 26.8, 2…
## $ Carbon_Sequestration_value <dbl> 2.39, 2.03, 0.96, 0.39, 2.08, 2.00, 1.74, …
## $ Stormwater_ft              <dbl> 150.3, 133.0, 35.4, 38.3, 115.5, 85.3, 99.…
## $ Stormwater_value           <dbl> 10.05, 8.89, 2.37, 2.56, 7.72, 5.70, 6.66,…
## $ Pollution_Removal_value    <dbl> 18.25, 16.15, 4.31, 4.65, 14.03, 10.36, 12…
## $ Pollution_Removal_oz       <dbl> 40.4, 35.8, 9.5, 10.3, 31.1, 22.9, 26.8, 2…
## $ Total_Annual_Benefits      <dbl> 30.69, 27.06, 7.64, 7.59, 23.82, 18.06, 20…
## $ Family                     <chr> "Pinaceae", "Pinaceae", "Rosaceae", "Fagac…
## $ Genus                      <chr> "Pseudotsuga", "Pseudotsuga", "Crataegus",…
## $ Common_name                <chr> "Douglas-fir", "Douglas-fir", "Lavalle haw…
## $ Genus_species              <chr> "Pseudotsuga menziesii", "Pseudotsuga menz…
## $ Functional_type            <chr> "CE", "CE", "BD", "BD", "CE", "CE", "CE", …
## $ Size                       <chr> "L", "L", "S", "L", "L", "L", "L", "L", "L…
## $ Native                     <chr> "Yes", "Yes", "No", "No", "Yes", "Yes", "Y…
## $ Edible                     <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ Nuisance                   <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ Origin                     <chr> "North America - from British Columbia sou…
## $ Species_factoid            <chr> "Bracts on cones look like a mouse's feet …

trees_raw %>% ggplot(mapping=aes(x=TreeHeight)) + geom_histogram(bins=50) + ggtitle("Histogram of park tree height")

## Warning: Removed 1 rows containing non-finite values (stat_bin).

Adding trees and river to the map

pdx_boundaries %>% ggplot() +
  geom_sf() +
  stat_density_2d(data=trees_raw, aes(X, Y))+
  geom_sf(data=river_boundaries, fill="blue", size=0.0)

# Converting coordinates for the park

p3 <- data.frame(lon=c(-122.675, -122.665), lat=c(45.48, 45.47))
sf_project(
  sf::st_crs(4326), # to: default datum used by geom_sf
  sf::st_crs(3857), # from: datum used by Portland Maps data
  p3
)

##           [,1]    [,2]
## [1,] -13656119 5697407
## [2,] -13655005 5695819

Mapping the trees in the park

pdx_boundaries %>% ggplot() +
  geom_sf(fill='palegoldenrod') +
  geom_sf(data=river_boundaries, fill="blue", size=0.0) +
  geom_sf(data=trees_shape) + 
  coord_sf(xlim=c( -13656119, -13655005), ylim=c(5695819, 5697407))+
  ggtitle("Wilamette Park")+
  theme_economist()

# Count of trees in the park

trees_in_park <- st_crop(trees_shape, xmin=-13656119, ymin=5695819, xmax=-13655005, ymax=5697407)

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

paste0("Total number of trees in original dataset: ", length(trees_shape$geometry))

## [1] "Total number of trees in original dataset: 25534"

paste0("Number of trees in our bounding box: ", length(trees_in_park$geometry))

## [1] "Number of trees in our bounding box: 370"

Challenge 6

Use OpenStreetMap tiles

Standard Map

park_base_map <- get_stamenmap(
  c(left=-122.674, bottom=45.47, right=-122.666, top=45.48), zoom=16, maptype="terrain")

## Source : http://tile.stamen.com/terrain/16/10435/23450.png

## Source : http://tile.stamen.com/terrain/16/10436/23450.png

## Source : http://tile.stamen.com/terrain/16/10437/23450.png

## Source : http://tile.stamen.com/terrain/16/10435/23451.png

## Source : http://tile.stamen.com/terrain/16/10436/23451.png

## Source : http://tile.stamen.com/terrain/16/10437/23451.png

## Source : http://tile.stamen.com/terrain/16/10435/23452.png

## Source : http://tile.stamen.com/terrain/16/10436/23452.png

## Source : http://tile.stamen.com/terrain/16/10437/23452.png

## Source : http://tile.stamen.com/terrain/16/10435/23453.png

## Source : http://tile.stamen.com/terrain/16/10436/23453.png

## Source : http://tile.stamen.com/terrain/16/10437/23453.png

ggmap(park_base_map)+
  ggtitle("Wilamette Park")

Watercolor map

park_base_map <- get_stamenmap(
  c(left=-122.674, bottom=45.47, right=-122.666, top=45.48), zoom=16, maptype="watercolor")

## Source : http://tile.stamen.com/watercolor/16/10435/23450.jpg

## Source : http://tile.stamen.com/watercolor/16/10436/23450.jpg

## Source : http://tile.stamen.com/watercolor/16/10437/23450.jpg

## Source : http://tile.stamen.com/watercolor/16/10435/23451.jpg

## Source : http://tile.stamen.com/watercolor/16/10436/23451.jpg

## Source : http://tile.stamen.com/watercolor/16/10437/23451.jpg

## Source : http://tile.stamen.com/watercolor/16/10435/23452.jpg

## Source : http://tile.stamen.com/watercolor/16/10436/23452.jpg

## Source : http://tile.stamen.com/watercolor/16/10437/23452.jpg

## Source : http://tile.stamen.com/watercolor/16/10435/23453.jpg

## Source : http://tile.stamen.com/watercolor/16/10436/23453.jpg

## Source : http://tile.stamen.com/watercolor/16/10437/23453.jpg

ggmap(park_base_map)+
  ggtitle("Wilamette Park")

# Adding trees

trees_park_projected <- st_transform(trees_in_park, sf::st_crs(4326))
ggmap(park_base_map) + 
  ggtitle("Wilamette Park")+
  geom_sf(data=trees_park_projected, inherit.aes = FALSE)

## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

#Coloring the trees and changing color option on the map

park_base_map <- get_stamenmap(
  c(left=-122.674, bottom=45.47, right=-122.666, top=45.48), zoom=16, maptype="terrain")

ggmap(park_base_map) + 
  ggtitle("Wilamette Park")+
  geom_sf(data=trees_park_projected, inherit.aes = FALSE,
          alpha=0.05,
          mapping=aes(color=Condition))

## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

# making a hybrid map

leaflet(trees_park_projected) %>% # we can pass in a dataframe to use...
  addEsriBasemapLayer(esriBasemapLayers$Imagery) %>%
  setView(-122.670, 45.475, zoom=16) %>%
  addCircleMarkers(label = ~Condition)

Challenge 7

A map using a different census column (or multiple columns!)

For this exercise I wanted to show density of people living alone in the household. Might be useful for a dating site.

pdx_lonesome <- read_excel("data/new_census.xlsx",
                              sheet="Census_2010_Neighborhoods", 
                              range="A7:D101", 
                              col_names = c("Neighborhood", "Total.Pop", "MaleAlone", "FemaleAlone")) %>%                    mutate(Neighborhood=as.factor(Neighborhood))

pdx_lonesome <- pdx_lonesome %>% 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"
))

boundaries_with_lonesome <- left_join(pdx_boundaries, pdx_lonesome, by=c("NAME"="Neighborhood"))

## Warning: Column `NAME`/`Neighborhood` joining factors with different levels,
## coercing to character vector

Bulding area map with high concentration of males living alone

boundaries_with_lonesome_males <- boundaries_with_lonesome %>% 
  mutate(adj_male_pop = (MaleAlone / (Shape_Area / 1E6))) %>% arrange(adj_male_pop)

boundaries_with_lonesome_males %>% ggplot() + 
  geom_sf(aes(fill=adj_male_pop)) +
  scale_fill_gradient() +
  geom_sf(data=river_boundaries, fill="turquoise3", size=0.0)+
  ggtitle("Neighborhoods with high concentration males living alone")

Building area map with high population of females living alone

boundaries_with_lonesome_females <-boundaries_with_lonesome %>% 
  mutate(adj_fem_pop = (FemaleAlone / (Shape_Area / 1E6))) %>% arrange(adj_fem_pop)

boundaries_with_lonesome_females %>% ggplot() + 
  geom_sf(aes(fill=adj_fem_pop)) +
  scale_fill_gradient() +
  geom_sf(data=river_boundaries, fill="turquoise3", size=0.0)+
  ggtitle("Neighborhoods with high concentration females living alone")

Challenge 8

Exploring color scales

boundaries_with_lonesome_females %>% ggplot(mapping=aes(x=adj_fem_pop))+
  geom_histogram(bins=50) + 
  ggtitle("Histogram of lonesome females")

## Warning: Removed 4 rows containing non-finite values (stat_bin).

min.pop <- min(boundaries_with_lonesome_females$adj_fem_pop, na.rm=TRUE)
max.pop <- max(boundaries_with_lonesome_females$adj_fem_pop, na.rm=TRUE)
diff.pop <- max.pop - min.pop
std.dev.pop <- sd(boundaries_with_lonesome_females$adj_fem_pop, na.rm = TRUE)

# some possible color break points
equal.interval <- seq(min.pop, max.pop, by= diff.pop / 4) 
quantile.interval <- quantile(boundaries_with_lonesome_females$adj_fem_pop, probs=seq(0, 1, by = 1/4), na.rm = TRUE)
std.interval <- c(seq(min.pop, max.pop, by=std.dev.pop), max.pop)
jenks.interval <- classIntervals(boundaries_with_lonesome_females$adj_fem_pop, n=6, style='jenks')$brks

## Warning in classIntervals(boundaries_with_lonesome_females$adj_fem_pop, : var
## has missing values, omitted in finding classes

boundaries_with_lonesome_females$pop.equal = cut(boundaries_with_lonesome_females$adj_fem_pop, breaks=equal.interval, include.lowest = TRUE)
boundaries_with_lonesome_females$pop.quantile = cut(boundaries_with_lonesome_females$adj_fem_pop, breaks=quantile.interval, include.lowest = TRUE)
boundaries_with_lonesome_females$pop.std = cut(boundaries_with_lonesome_females$adj_fem_pop, breaks=std.interval, include.lowest = TRUE)
boundaries_with_lonesome_females$pop.jenks = cut(boundaries_with_lonesome_females$adj_fem_pop, breaks=jenks.interval, include.lowest = TRUE)

popDensBarChart <- function(break_col) {
  boundaries_with_lonesome_females %>% 
  filter(!is.na(adj_fem_pop)) %>% 
  ggplot(mapping = aes(x=fct_reorder(MAPLABEL, -adj_fem_pop), y=adj_fem_pop)) +
    geom_col(aes(fill=.data[[break_col]])) +
    # theme(axis.text.x=element_text(angle=270, hjust=0, size=6, vjust=0.5)) +
  theme(axis.text.x=element_blank(), axis.ticks.x=element_blank()) +
    scale_y_continuous(name=expression("People/km" ^2)) +
  scale_x_discrete(name=NULL) +
    scale_fill_discrete(guide=FALSE) +
  ggtitle(break_col)
}

It seems like using 4

plot_grid(
  popDensBarChart("pop.equal"),
  popDensBarChart("pop.quantile"),
  popDensBarChart("pop.std"),
  popDensBarChart("pop.jenks"),
  nrow=2, ncol=2
)

# Final graphs

pop.plot <- function(col){
  boundaries_with_lonesome_females %>%
  ggplot()+
  geom_sf(aes(fill=.data[[col]]))+
  scale_fill_brewer(palette="YlOrBr", name=expression("people/km" ^2))+
  geom_sf(data=river_boundaries, fill='turquoise3', size=0.0)+
  ggtitle(col)
}

pop.plot('pop.equal')

  pop.plot("pop.quantile")

  pop.plot("pop.std")

  pop.plot("pop.jenks")