Preliminary tree_year and demo_canopy maps

scripts/demographics_canopy_map.R

@@ -0,0 +1,91 @@
+# Title     : demographics_canopy_maps.R
+# Objective : compare demographics at census tract level with tree equity
+# Created by: Charlotte
+# Created on: 7 Feb 2021
+
+#Remaining questions:
+# Total population of different categories aren't the same? which to use?
+# Remove census areas that are over cemetery and airport?
+
+#TO DO:
+#Get neighborhood names for tract instead of geo_id
+#make more generalizable and less hard coded
+#add canopy coverage as layer (got error when tried to, maybe too large)
+#Check which libraries are needed
+
+library(tidyverse)
+library(leaflet)
+library(RColorBrewer)
+library(htmltools)
+library(htmlwidgets)
+library(sf)
+library(raster)
+library("rgdal")
+library(units)
+library(lwgeom)
+
+#get data and functions
+#setwd('Path_to/EcoAction')
+source("src/load_data.R")
+source("src/demos_canopy_utils.R")
+#separately, map demography per area or population and canopy cover per area
+
+#For tracts-get demographics, canopy, and make map
+#get demographics
+demo_tract <- read_demographics_csv('tract')
+#get canopy area
+canopy_tract <- read_land_area_csv('tract')
+#read geos of tracts and add to demo_tract and canopy_tract
+tracts <- read_geos_tract()
+#add geometry
+demo_tract <- merge(demo_tract, tracts, by="geo_id")
+canopy_tract <- merge(canopy_tract, tracts, by="geo_id")
+
+#add population density to demo_tract
+demo_tract$pop_density = demo_tract$tot_pop_race/canopy_tract$area_m_sq
+
+#equity measures: canopy area by person and canopy area by person in poverty
+demo_tract$can_pop = canopy_tract$area_canopy/demo_tract$tot_pop_race
+demo_tract$can_pov = canopy_tract$area_canopy/demo_tract$pop_in_poverty
+
+#make maps with different demographic and canopy values for each
+#tract (can overlay layers in secon map)
+tract_demo_canopy_base <- make_base_groups(demo_tract, canopy_tract, 'tract_demo_canopy_base.html')
+tract_demo_canopy_overlay <- make_overlay_groups(demo_tract, canopy_tract, 'tract_demo_canopy_overlay.html')
+
+#pairwise plots
+corr_df = data.frame(demo_tract$pct_in_poverty, demo_tract$pct_nonwhite,
+               canopy_tract$pct_plantable, canopy_tract$pct_open_plantable,
+               canopy_tract$pct_canopy)
+pairs(corr_df, pch=18)
+
+#make same maps at block_group level (smaller than tracts)
+#get demographics
+demo_bg <- read_demographics_csv('block_group')
+#get canopy area
+canopy_bg <- read_land_area_csv('block_group')
+#read geos of bgs and add to demo_bg and canopy_bg
+block_group <- read_geos_block_group()
+#merge bgs with geometries
+demo_bg <- merge(demo_bg, block_group, by="geo_id")
+canopy_bg <- merge(canopy_bg, block_group, by="geo_id")
+
+#add population density to demo_bg
+demo_bg$pop_density = demo_bg$tot_pop_race/canopy_bg$area_m_sq
+
+#equity measures: canopy area by person and canopy area by person in poverty
+demo_bg$can_pop = canopy_bg$area_canopy/demo_bg$tot_pop_race
+demo_bg$can_pov = canopy_bg$area_canopy/demo_bg$pop_in_poverty
+
+#make maps with different demographic and canopy values for each
+#block group (can overlay layers in secon map)
+bg_demo_canopy_base <- make_base_groups(demo_bg, canopy_bg, 'bg_demo_canopy_base.html')
+bg_demo_canopy_overlay <- make_overlay_groups(demo_bg, canopy_bg, 'bg_demo_canopy_overlay.html')
+
+#pairwise graphs comparing some different values in each block group
+corr_df = data.frame(demo_bg$pct_in_poverty, demo_bg$pct_nonwhite,
+                     canopy_bg$pct_plantable, canopy_bg$pct_open_plantable,
+                     canopy_bg$pct_canopy)
+pairs(corr_df, pch=18)
+
+plot(demo_bg$pct_in_poverty, canopy_bg$pct_canopy, pch=18)

scripts/map_year.R

@@ -0,0 +1,55 @@
+# Title     : map_year.R
+# Objective : Make maps of trees planted by Eco Arlington per year
+# Created by: Charlotte
+# Created on: 7 Feb 2021
+
+#TO DO: make some tree planted per census area map?
+
+library(tidyverse)
+library(leaflet)
+library(RColorBrewer)
+library(htmltools)
+library(htmlwidgets)
+library(sf)
+library(raster)
+library("rgdal")
+library(units)
+library(lwgeom)
+library(ggplot2)
+
+
+#get data and functions
+#setwd('Path_to/EcoAction')
+source("src/load_data.R")
+source("src/map_year_utils.R")
+trees_sf <- read_tree_data_subset()
+#if you get "'data/tree_data_consolidated - trees.csv' does not exist in
+#current working directory" error, ensure you've gotten file from google drive
+
+#for better labeling add jitter to points. This means it is possible to see
+#labels of trees at the same exact location.
+#note that this means the locations of trees is no longer exact
+trees_sf_jitter <- trees_sf
+trees_sf_jitter$geometry <- st_jitter(trees_sf_jitter$geometry,
+                                      amount = 0.0001, factor = 0.0001)
+
+#make map without jitter
+ym <- year_map(trees_sf, "map_trees_by_year.html")
+ylm <- year_layer_map(trees_sf, "map_trees_layer_by_year.html")
+
+#make map with jitter
+ym_jitter <- year_map(trees_sf_jitter, "map_trees_jitter_by_year.html")
+ylm_jitter <- year_layer_map(trees_sf_jitter, "map_trees_jitter_layer_by_year.html")
+
+
+#histogram plots to better understand year and tree type patterns
+#number of trees planted each year
+y <-hist(trees_sf$year)
+#number of trees planted by type of tree, color-coded by year planted
+year_type <-trees_sf%>%
+    group_by(year)%>%
+    count(tree_name)
+year_type_plot <- ggplot(year_type, aes(fill=year, y=n, x=tree_name)) +
+    geom_bar(position="stack", stat="identity")
+year_type_plot + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
+

src/demos_canopy_utils.R

@@ -0,0 +1,189 @@
+# Title     : demographics_canopy_utils.R
+# Objective : functions to make maps displaying various demographic and
+#             canopy values across census areas (tract or block group)
+# Created by: Charlotte
+# Created on: 15 Feb 2021
+
+#TO Do:
+#Add legends?
+
+library(tidyverse)
+library(leaflet)
+library(RColorBrewer)
+library(htmltools)
+library(htmlwidgets)
+library(sf)
+library(raster)
+library("rgdal")
+library(units)
+library(lwgeom)
+
+#helper function to make labels for each census area (hard-coded value types)
+make_labels <- function(demo, canopy){
+  labels <- sprintf(
+    "<strong>Geo ID: %s</strong>
+    <br/>%g percent canopy cover
+    <br/>%g percent open plantable land
+    <br/>%g total population
+    <br/>%g population density (people per square meter)
+    <br/>%g percent Black, Asian, Pacific Islander, American Indian,
+    <br/> Other, or two or more races
+    <br/>%g percent population below poverty level
+    <br/>%g canopy cover per person (square meter per person)
+    <br/>%g canopy cover per person below poverty
+    <br/>level (square meter per person below  poverty level)",
+    canopy$geo_id, canopy$pct_canopy,
+    canopy$pct_open_plantable, demo$tot_pop_race,
+    demo$pop_density,
+    demo$pct_nonwhite, demo$pct_in_poverty,
+    demo$can_pop, demo$can_pov
+  ) %>% lapply(htmltools::HTML)
+
+  return(labels)
+}
+
+# helper function to make basic map with hard-coded demographic and canopy
+# value types
+make_map_basic <- function(demo, canopy,labels, layers){
+
+
+    #make color ramp for each demo/canopy value type
+    qpal <- leaflet::colorQuantile(
+        palette = "RdPu",
+        domain = demo$pct_nonwhite, n=7)
+    qpal_can <- leaflet::colorQuantile(
+        palette = "YlGn",
+        domain = canopy$pct_canopy, n=7)
+    qpal_openplant <- leaflet::colorQuantile(
+        palette = "GnBu",
+        domain = canopy$pct_open_plantable, n=7)
+    qpal_pov <- leaflet::colorQuantile(
+        palette = "YlOrRd",
+        domain = demo$pct_in_poverty, n=7)
+    qpal_can_pov <- leaflet::colorQuantile(
+        palette = "YlOrRd",
+        domain = demo$can_pov, n=7)
+    qpal_density <- leaflet::colorQuantile(
+        palette = "Greys",
+        domain = demo$pop_density, n=7)
+    qpal_can_pop <- leaflet::colorQuantile(
+        palette = "Greys",
+        domain = demo$can_pop, n=7)
+
+    #make map with colored polygons for each value type + census area
+    map <- demo$geometry %>% leaflet() %>%
+        # Base groups
+        addProviderTiles("CartoDB.Positron")%>%
+        #addProviderTiles("CartoDB.Positron", group = 'blank')%>%
+        # Overlay groups
+        addPolygons(color = qpal_can(canopy$pct_canopy) ,
+                    weight=1,
+                    smoothFactor=0.5,
+                    opacity=1.0,
+                    fillOpacity=0.5,
+                    group = layers[1],
+                    label = labels
+        )%>%
+        addPolygons(color = qpal_openplant(canopy$pct_open_plantable) ,
+                    weight=1,
+                    smoothFactor=0.5,
+                    opacity=1.0,
+                    fillOpacity=0.5,
+                    group = layers[2],
+                    label = labels
+        )%>%
+        addPolygons(color = qpal_density(demo$pop_density) ,
+                    weight=1,
+                    smoothFactor=0.5,
+                    opacity=1.0,
+                    fillOpacity=0.5,
+                    group = layers[3],
+                    label = labels
+        )%>%
+        addPolygons(color = qpal(demo$pct_nonwhite) ,
+                    weight=1,
+                    smoothFactor=0.5,
+                    opacity=1.0,
+                    fillOpacity=0.5,
+                    group = layers[4],
+                    label = labels
+        )%>%
+        addPolygons(color = qpal_pov(demo$pct_in_poverty) ,
+                    weight=1,
+                    smoothFactor=0.5,
+                    opacity=1.0,
+                    fillOpacity=0.5,
+                    group = layers[5],
+                    label = labels
+        )%>%
+        addPolygons(color = qpal_can_pop(demo$can_pop) ,
+                    weight=1,
+                    smoothFactor=0.5,
+                    opacity=1.0,
+                    fillOpacity=0.5,
+                    group = layers[6],
+                    label = labels
+        )%>%
+        addPolygons(color = qpal_can_pov(demo$can_pov) ,
+                    weight=1,
+                    smoothFactor=0.5,
+                    opacity=1.0,
+                    fillOpacity=0.5,
+                    group = layers[7],
+                    label = labels
+        )
+    return(map)
+}
+
+# with this map type, only one layer (demo/canopy value type) is displayed
+# at one time
+make_base_groups <- function(demo, canopy, map_name){
+
+  #make layers
+    layers = c("percent canopy cover",
+               "percent open plantable land",
+               "population density (people per square meter)",
+               "percent Black, Asian, Pacific Islander, American Indian, <br/> Other, or two or more races",
+               "percent population below poverty level",
+               "canopy cover per person (square meter per person)",
+               "canopy cover per person below poverty <br/>level (square meter per person below  poverty level)")
+  #make labels for each census area
+  labels <- make_labels(demo, canopy)
+  #make basic map
+  map <- make_map_basic(demo, canopy, labels, layers)
+  #add layer control to map
+  map <- map%>%
+    addLayersControl(
+      baseGroups = layers,
+      options = layersControlOptions(collapsed = FALSE)
+    )
+  # save to file
+  htmlwidgets::saveWidget(map, file = map_name, selfcontained=TRUE)
+    return(map)
+}
+
+# with this map type, several layers (demo/canopy value type) can be layered
+# together
+make_overlay_groups <- function(demo, canopy, map_name){
+    #make layers
+    layers = c("percent canopy cover",
+               "percent open plantable land",
+               "population density (people per square meter)",
+               "percent Black, Asian, Pacific Islander, American Indian, <br/> Other, or two or more races",
+               "percent population below poverty level",
+               "canopy cover per person (square meter per person)",
+               "canopy cover per person below poverty <br/>level (square meter per person below  poverty level)")
+    #make labels for each census area
+    labels <- make_labels(demo, canopy)
+    #make basic map
+    map <- make_map_basic(demo, canopy, labels, layers)
+    #add layer control to map
+    map <- map%>%
+        addLayersControl(
+            overlayGroups = layers,
+            options = layersControlOptions(collapsed = FALSE)
+        )
+    # save to file
+    htmlwidgets::saveWidget(map, file = map_name, selfcontained=TRUE)
+    return(map)
+}