-
Notifications
You must be signed in to change notification settings - Fork 0
/
boulder_trees.Rmd
123 lines (106 loc) · 6.42 KB
/
boulder_trees.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
---
title: "Boulder Tree Canopy"
author: "Jacob Paul and Jack Sandberg"
date: "10/3/2018"
output:
html_document: default
pdf_document: default
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(dplyr)
library(leaflet)
library(rgdal)
library(ggplot2)
library(rgeos)
library(raster)
```
###Import the Data
```{r}
shape <- readOGR('./TreesCityBoulder', layer='TreesCityBoulder')
shape_latlon <- spTransform(shape, CRS("+proj=longlat +datum=WGS84"))
```
#Preliminary Exploration
##Looking at Groups of the EAB column
We're not sure what exactly these groupings are and what they represent, but feel that they might be useful to some of our analysis.
Below, we mapped out the groups to see if we could find a pattern, and the assigned Group seems to depend on the region of the city that the tree is in.
```{r}
filtered <- shape_latlon[shape_latlon$EAB != 'Unassigned',]
pal = colorFactor(palette = "RdYlBu", domain=filtered$EAB)
leaflet(data=filtered) %>% addProviderTiles(providers$CartoDB.Positron) %>%
addCircles(color=~pal(filtered$EAB)) %>% addLegend("bottomright", pal=pal, filtered$EAB)
```
##Finding the Ash Trees
Below is a map of all of the ash trees documented by the city labled by the kind of ash tree (white or green). We determined from this map that there is no evident pattern of ash tree existence just by looking at this map.
```{r}
filtered <- shape_latlon[shape_latlon$COMMONNAME == "Green Ash" | shape_latlon$COMMONNAME == "White Ash",]
pal = colorFactor(palette = "RdYlBu", domain=filtered$COMMONNAME)
leaflet(data=filtered) %>% addProviderTiles(providers$CartoDB.Positron) %>%
addCircles(color=~pal(filtered$COMMONNAME)) %>% addLegend("bottomright", pal=pal, filtered$COMMONNAME)
```
##Mapping Tree Conditions
Although the color scheme is off, we were able to determine that there are no surface level patterns that can be drawn from this map either.
```{r}
filtered <- shape_latlon[shape_latlon$COMMONNAME == "Green Ash" | shape_latlon$COMMONNAME == "White Ash",]
pal = colorFactor(palette = "RdYlBu", domain=c("Dead","Very Poor", "Poor", "Fair", "Good", "Excellent", "N/A"))
leaflet(data=filtered) %>% addProviderTiles(providers$CartoDB.Positron) %>%
addCircles(color=~pal(filtered$CONDITION)) %>% addLegend("bottomright", pal=pal, filtered$CONDITION)
```
##Diving into Ash Tree grouping
The heatmap below shows areas where the relative amount of ash trees is high. It's clear that some parts of the city contain a lot of ash trees, but we still don't know how these numbers compare to the total number of trees.
```{r}
library(leaflet.extras)
shape_latlon$is_ash <- shape_latlon$COMMONNAME == "Green Ash" | shape_latlon$COMMONNAME == "White Ash"
pal = colorFactor(palette = topo.colors(2), domain=c("TRUE", "FALSE"))
leaflet(data=filtered) %>% addProviderTiles(providers$CartoDB.Positron) %>% addHeatmap(radius=6)
#addCircles(color=~pal(shape_latlon$is_ash)) %>%
#addLegend("bottomright", pal=pal, shape_latlon$is_ash)
```
#Grouping by census block
We needed a way to group trees into different segments. Ideally, we would group trees into a grid format because census blocks are a way for grouping people rather than trees, but this will work for now. First, we need to import the census blocks.
```{r}
census_blocks <- readOGR('./ACS1216_bg', layer='ACS1216_bg')
census_blocks_latlong <- spTransform(census_blocks, CRS("+proj=longlat +datum=WGS84"))
```
The shapefile above contains census blocks for the entire state of Colorado, and we only need the census blocks in Boulder. We used the Boulder City Limits shapefile avaialable on the City of Boulder's Open Data Catalog.
```{r}
citylims <- readOGR('./Boulder_City_Limits', layer='Boulder_City_Limits')
citylims_latlong <- spTransform(citylims, CRS("+proj=longlat +datum=WGS84"))
#leaflet(data=citylims_latlong) %>% addProviderTiles(providers$CartoDB.Positron) %>% addPolygons()
```
Now, we can use the city limits shape file to crop our census blocks down to just those within the city limits.
```{r}
filtered_blocks <- crop(census_blocks_latlong, citylims_latlong)
blocks_objectid <- filtered_blocks[c("OBJECTID")]
```
The OBJECTID column in the census blocks data will serve as a good value to group our trees by.
```{r}
block_ref <- intersect(shape_latlon, blocks_objectid)
summary(block_ref$d)
```
These lines of code were used to verify that the grouping occured succesfully. They have been commented out now becuase the map was unecessary to keep in the analysis.
```{r}
#pal = colorFactor(palette = "RdYlBu", domain=block_ref$d)
#leaflet(data=block_ref[block_ref$COMMONNAME == "White Ash" | block_ref$COMMONNAME == "Green Ash",]) %>% #addProviderTiles(providers$CartoDB.Positron) %>% addPolygons(data=filtered_blocks) %>% addCircles(data = block_ref, color = #~pal(block_ref$d))
```
The goal is to use z-score to determine the percentage of ash trees in each block relative to the average. To make sure that this is a valid method of analysis, we plotted the percentage of ash trees in each block to make sure they roughly followed a normal distribution. Since the data follow a normal distribution, we calculated the z-score for each block.
```{r}
#data.frame(block_ref)
block_gb <- data.frame(block_ref) %>% group_by(d) %>% tally()
block_gb_ash <- data.frame(block_ref[block_ref$COMMONNAME == "White Ash" | block_ref$COMMONNAME == "Green Ash",]) %>% group_by(d) %>% tally()
block_gb_merge <- merge(block_gb, block_gb_ash, by="d")
block_gb_merge$ratio <- block_gb_merge$n.y/block_gb_merge$n.x
qplot(block_gb_merge$ratio, geom="histogram", bins=12)
block_gb_merge$zscore <- (block_gb_merge$ratio-mean(block_gb_merge$ratio))/sd(block_gb_merge$ratio)
```
In order to get a better idea what this meant spatially, we plotted these blocks on a map, color coded according to their z-score. Here we can more easily spot where there are statistically more ash trees than average. This can be used to focus EAB mitigation efforts as well as to anticipate where the largest impact will be.
```{r}
sp_blocks_zscore <- merge(blocks_objectid, block_gb_merge, by.x="OBJECTID", by.y="d")
colors <- colorBin("PRGn", domain = sp_blocks_zscore$zscore, bins = c(-3,-2,-1,0,1,2,3))
leaflet(sp_blocks_zscore) %>% addProviderTiles(providers$CartoDB.Positron) %>% addPolygons(color=colors(sp_blocks_zscore$zscore), weight=3, fillOpacity = .7) %>% addLegend("bottomright",
pal = colors,
values = sp_blocks_zscore$zscore,
title = "Z-score of % of ash trees")# %>% addHeatmap(data=filtered, radius=6)
```
```{r}
```