1
+
2
+ library(plotly )
3
+ library(leaflet )
4
+ library(plyr )
5
+ library(dplyr )
6
+ library(ggplot2 )
7
+ library(rgdal )
8
+ library(tidyverse )
9
+ library(sf )
10
+ library(rgeos )
11
+ library(httr )
12
+ library(jsonlite )
13
+ library(raster )
14
+ library(data.table )
15
+ library(leaflet.providers )
16
+ library(knitr )
17
+ library(kableExtra )
18
+ library(reactable )
19
+ library(readxl )
20
+ library(formattable )
21
+ library(reactablefmtr )
22
+
23
+
24
+ col_BoldRiverBlue = " #242456"
25
+ col_BoldSunYellow = " #FFD450"
26
+ col_BoldGrassGreen = " #2A553E"
27
+ col_BoldEarthGrey = " #7B7A66"
28
+ col_BoldBrickOrange = " #F26640"
29
+ col_LightRiverBlue = " #E3E6FF"
30
+ col_LightSunYellow = " #FFEDBA"
31
+ col_LightGrassGreen = " #C4F4D5"
32
+ col_LightEarthGrey = " #ECE2D6"
33
+ col_LightBrickOrange = " #FED3CF"
34
+ col_White = " #FFFFFF"
35
+ col_Black = " #000000"
36
+
37
+
38
+ # read Boundary
39
+ boundary = st_read(" https://cities-urbanshift.s3.eu-west-3.amazonaws.com/data/boundaries/CRI-San_Jose-boundary.geojson" , quiet = TRUE )
40
+
41
+ boundary_municipality = st_read(" ./data/biodiversity/boundaries-CRI-San_Jose-ADM2.geojson" ,
42
+ quiet = TRUE )
43
+
44
+
45
+ biodiversity_baseline_indicators <- read_excel(" ./data/biodiversity baseline indicators.xlsx" , sheet = " SCORES" )
46
+
47
+ boundary_municipality $ Municipality = str_remove(boundary_municipality $ shapeName.1 , " Cantón " )
48
+ # change encoding
49
+ boundary_municipality $ Municipality = iconv(boundary_municipality $ Municipality ,from = " UTF-8" ,to = " ASCII//TRANSLIT" )
50
+
51
+ # join with geo
52
+ biodiversity_baseline_indicators_geo = boundary_municipality %> %
53
+ left_join(biodiversity_baseline_indicators ,
54
+ by = c(" Municipality" )) %> %
55
+ rename(I1_percent_natural_areas_value = ' I-1 raw' ,
56
+ I1_percent_natural_areas_score = ' I-1 score' ,
57
+ I2_connectivity_value = ' I-2 raw' ,
58
+ I2_connectivity_score = ' I-2 score' ,
59
+ I3_percent_birds_built_area_value = ' I-3 raw' ,
60
+ I3_percent_birds_built_area_score = ' I-3 score' ,
61
+ I8_percent_protected_area_value = ' I-8 raw' ,
62
+ I8_percent_protected_area_score = ' I-8 score' ,
63
+ I10_water_value = ' I-10 raw' ,
64
+ I10_water_score = ' I-10 score' ,
65
+ I12_recreational_services_value = ' I-12 raw' ,
66
+ I12_recreational_services_score = ' I-12 score' ,
67
+ I13_proximity_parks_value = ' I-13 raw' ,
68
+ I13_proximity_parks_score = ' I-13 score' ) %> %
69
+ mutate_at(" I1_percent_natural_areas_score" , as.character ) %> %
70
+ mutate_at(" I2_connectivity_score" , as.character ) %> %
71
+ mutate_at(" I3_percent_birds_built_area_score" , as.character )%> %
72
+ mutate_at(" I8_percent_protected_area_score" , as.character ) %> %
73
+ mutate_at(" I10_water_score" , as.character ) %> %
74
+ mutate_at(" I12_recreational_services_score" , as.character ) %> %
75
+ mutate_at(" I13_proximity_parks_score" , as.character )
76
+
77
+ # plot matrix
78
+
79
+
80
+
81
+ biodiversity_baseline_scores = biodiversity_baseline_indicators_geo %> %
82
+ as.data.frame() %> %
83
+ drop_na(I1_percent_natural_areas_value ) %> %
84
+ mutate(I1_percent_natural_areas_score = as.numeric(I1_percent_natural_areas_score ),
85
+ I2_connectivity_score = as.numeric(I2_connectivity_score ),
86
+ I3_percent_birds_built_area_score = as.numeric(I3_percent_birds_built_area_score ),
87
+ I8_percent_protected_area_score = as.numeric(I8_percent_protected_area_score ),
88
+ I10_water_score = as.numeric(I10_water_score ),
89
+ I12_recreational_services_score = as.numeric(I12_recreational_services_score ),
90
+ I13_proximity_parks_score = as.numeric(I13_proximity_parks_score )) %> %
91
+ dplyr :: select(Municipality ,
92
+ " I1 - Percent of natural areas" = I1_percent_natural_areas_score ,
93
+ " I2 - Connectivity" = I2_connectivity_score ,
94
+ " I3 - Birds in built areas" = I3_percent_birds_built_area_score ,
95
+ " I8 - Protected areas" = I8_percent_protected_area_score ,
96
+ " I10 - Water" = I10_water_score ,
97
+ " I12 - Recreational services" = I12_recreational_services_score ,
98
+ " I13 - Proximity to parks" = I13_proximity_parks_score ) %> %
99
+ mutate(" Biodiversity Index" = rowSums(. [2 : 8 ])) %> %
100
+ mutate_at(vars(" I1 - Percent of natural areas" : " I13 - Proximity to parks" ), ~ cell_spec(
101
+ . , " html" ,
102
+ background = ifelse(. > = 4 , " green" , ifelse(. > = 3 , " yellowgreen" , ifelse(. > = 2 , " yellow" , ifelse(. > = 1 , " orange" , " red" ))))
103
+ )) %> %
104
+ mutate_at(vars(" Biodiversity Index" ), ~ cell_spec(
105
+ . , " html" ,
106
+ font_size = " x-large" ,
107
+ background = ifelse(. > = 28 , " green" , ifelse(. > = 21 , " yellowgreen" , ifelse(. > = 14 , " yellow" , ifelse(. > = 7 , " orange" , " red" ))))
108
+ ))
109
+
110
+ biodiversity_baseline_scores %> %
111
+ kable(format = " html" , escape = FALSE ) %> %
112
+ kable_styling(" striped" , full_width = FALSE ) %> %
113
+ scroll_box(width = " 100%" , height = " 700px" )
114
+
115
+
116
+
117
+ # map
118
+
119
+ biodiversity_baseline_scores = biodiversity_baseline_indicators_geo %> %
120
+ as.data.frame() %> %
121
+ drop_na(I1_percent_natural_areas_value ) %> %
122
+ mutate(I1_percent_natural_areas_score = as.numeric(I1_percent_natural_areas_score ),
123
+ I2_connectivity_score = as.numeric(I2_connectivity_score ),
124
+ I3_percent_birds_built_area_score = as.numeric(I3_percent_birds_built_area_score ),
125
+ I8_percent_protected_area_score = as.numeric(I8_percent_protected_area_score ),
126
+ I10_water_score = as.numeric(I10_water_score ),
127
+ I12_recreational_services_score = as.numeric(I12_recreational_services_score ),
128
+ I13_proximity_parks_score = as.numeric(I13_proximity_parks_score )) %> %
129
+ dplyr :: select(Municipality ,
130
+ " I1 - Percent of natural areas" = I1_percent_natural_areas_score ,
131
+ " I2 - Connectivity" = I2_connectivity_score ,
132
+ " I3 - Birds in built areas" = I3_percent_birds_built_area_score ,
133
+ " I8 - Protected areas" = I8_percent_protected_area_score ,
134
+ " I10 - Water" = I10_water_score ,
135
+ " I12 - Recreational services" = I12_recreational_services_score ,
136
+ " I13 - Proximity to parks" = I13_proximity_parks_score ) %> %
137
+ mutate(" Biodiversity Index" = rowSums(. [2 : 8 ]))
138
+
139
+ # join with geo
140
+ biodiversity_baseline_scores_geo = boundary_municipality %> %
141
+ left_join(biodiversity_baseline_scores ,
142
+ by = c(" Municipality" ))
143
+
144
+
145
+ # define color palette for I1 levels
146
+ pal_Index <- colorNumeric(palette = " Greens" ,
147
+ domain = biodiversity_baseline_scores_geo $ `Biodiversity Index` ,
148
+ na.color = " transparent" ,
149
+ revers = FALSE )
150
+
151
+ # define labels
152
+
153
+ labels_Index <- sprintf(" <strong>%s</strong><br/>%s: %s" ,
154
+ biodiversity_baseline_scores_geo $ shapeName.1 ,
155
+ " Biodiversity index" ,
156
+ biodiversity_baseline_scores_geo $ `Biodiversity Index` ) %> %
157
+ lapply(htmltools :: HTML )
158
+
159
+ # plot map
160
+ leaflet(height = 500 , width = " 100%" ) %> %
161
+ addTiles() %> %
162
+ addProviderTiles(" OpenStreetMap.France" , group = " OSM" ) %> %
163
+ addProviderTiles(providers $ CartoDB.DarkMatter , group = " CartoDB" ) %> %
164
+ addPolygons(data = biodiversity_baseline_scores_geo ,
165
+ group = " Biodiversity index" ,
166
+ fillColor = ~ pal_Index(`Biodiversity Index` ),
167
+ weight = 1 ,
168
+ opacity = 1 ,
169
+ color = " grey" ,
170
+ fillOpacity = 0.7 ,
171
+ label = labels_Index ,
172
+ highlightOptions = highlightOptions(color = " black" , weight = 2 ,
173
+ bringToFront = FALSE ),
174
+ labelOptions = labelOptions(
175
+ style = list (" font-weight" = " normal" , padding = " 3px 6px" ),
176
+ textsize = " 15px" ,
177
+ direction = " auto" )) %> %
178
+ addLegend(pal = pal_Index ,
179
+ values = biodiversity_baseline_scores_geo $ `Biodiversity Index` ,
180
+ opacity = 0.9 ,
181
+ title = " Biodiversity index (2020)" ,
182
+ position = " topright" ,
183
+ labFormat = labelFormat(suffix = " " )) %> %
184
+ # Layers control
185
+ addLayersControl(
186
+ baseGroups = c(" OSM" ," CartoDB" ),
187
+ overlayGroups = c(" Biodiversity index" ),
188
+ options = layersControlOptions(collapsed = FALSE )
189
+ )
190
+
191
+ # plot map -----
192
+
193
+ # define color palette for I1 levels
194
+ pal_I1 <- colorNumeric(palette = " Greens" ,
195
+ domain = biodiversity_baseline_indicators_geo $ I1_percent_natural_areas_value ,
196
+ na.color = " transparent" ,
197
+ revers = FALSE )
198
+
199
+ # define color palette for S1 levels
200
+ pal_score <- colorFactor(palette = c(" green" ," yellowgreen" ," yellow" ," orange" ," red" ),
201
+ levels = c(" 0" ," 1" ," 2" ," 3" ," 4" ),
202
+ na.color = " transparent" ,
203
+ revers = TRUE )
204
+
205
+ # define labels
206
+
207
+ labels_I1 <- sprintf(" <strong>%s</strong><br/>%s: %s %s<br/>%s: %s" ,
208
+ biodiversity_baseline_indicators_geo $ shapeName.1 ,
209
+ " Proportion of natural areas" ,
210
+ round(biodiversity_baseline_indicators_geo $ I1_percent_natural_areas_value , 2 ), " " ,
211
+ " Score" ,biodiversity_baseline_indicators_geo $ I1_percent_natural_areas_score ) %> %
212
+ lapply(htmltools :: HTML )
213
+
214
+
215
+
216
+ # plot map
217
+ leaflet(height = 500 , width = " 100%" ) %> %
218
+ addTiles() %> %
219
+ addProviderTiles(providers $ CartoDB.DarkMatter , group = " CartoDB" ) %> %
220
+ addProviderTiles(" OpenStreetMap.France" , group = " OSM" ) %> %
221
+ addPolygons(data = biodiversity_baseline_indicators_geo ,
222
+ group = " Porportion of natural areas" ,
223
+ fillColor = ~ pal_I1(I1_percent_natural_areas_value ),
224
+ weight = 1 ,
225
+ opacity = 1 ,
226
+ color = " grey" ,
227
+ fillOpacity = 0.7 ,
228
+ label = labels_I1 ,
229
+ highlightOptions = highlightOptions(color = " black" , weight = 2 ,
230
+ bringToFront = FALSE ),
231
+ labelOptions = labelOptions(
232
+ style = list (" font-weight" = " normal" , padding = " 3px 6px" ),
233
+ textsize = " 15px" ,
234
+ direction = " auto" )) %> %
235
+ addLegend(pal = pal_I1 ,
236
+ values = biodiversity_baseline_indicators_geo $ I1_percent_natural_areas_value ,
237
+ opacity = 0.9 ,
238
+ title = " % natural areas (2020)" ,
239
+ position = " topright" ,
240
+ labFormat = labelFormat(suffix = " " )) %> %
241
+ # I1 score layer
242
+ addPolygons(data = biodiversity_baseline_indicators_geo ,
243
+ group = " Porportion of natural areas Score" ,
244
+ fillColor = ~ pal_score(I1_percent_natural_areas_score ),
245
+ weight = 1 ,
246
+ opacity = 1 ,
247
+ color = " grey" ,
248
+ fillOpacity = 0.7 ,
249
+ label = labels_I1 ,
250
+ highlightOptions = highlightOptions(color = " black" , weight = 2 ,
251
+ bringToFront = FALSE ),
252
+ labelOptions = labelOptions(
253
+ style = list (" font-weight" = " normal" , padding = " 3px 6px" ),
254
+ textsize = " 15px" ,
255
+ direction = " auto" )) %> %
256
+ # I1 score legend
257
+ addLegend(colors = c(" green" ," yellowgreen" ," yellow" ," orange" ," red" ),
258
+ labels = c(" 4 (> 20.0%)" ,
259
+ " 3 (14.0% – 20.0%)" ,
260
+ " 2 (7.0% – 13.9%)" ,
261
+ " 1 (0% – 6.9%)" ,
262
+ " 0 (< 1.0%)" ),
263
+ opacity = 0.9 ,
264
+ title = " Natural areas scores (2020)" ,
265
+ position = " bottomleft" ,
266
+ labFormat = labelFormat(suffix = " " )) %> %
267
+ # Layers control
268
+ addLayersControl(
269
+ baseGroups = c(" CartoDB" ," OSM" ),
270
+ overlayGroups = c(" Porportion of natural areas" ," Porportion of natural areas Score" ),
271
+ options = layersControlOptions(collapsed = FALSE )
272
+ ) %> %
273
+ hideGroup(" Porportion of natural areas Score" )
0 commit comments