-
Notifications
You must be signed in to change notification settings - Fork 0
/
LST_MassachusettsPDF.Rmd
5155 lines (4646 loc) · 268 KB
/
LST_MassachusettsPDF.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
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
---
title: "Heat Risk in Massachusetts"
author: "Marcos Luna and Neenah Estrella-Luna"
date: "`r Sys.Date()`"
output:
bookdown::pdf_document2:
toc: true
toc_depth: 3
lot: yes
lof: yes
number_sections: true
fig_caption: yes
includes:
in_header: my_header.tex
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo=FALSE, message=FALSE, warning=FALSE)
```
\pagebreak
# Heat risk in Massachusetts
Heat, or hot weather, is the leading weather-related cause of death in the U.S.[^NWShazards] Between 1999 and 2010, the U.S. Centers for Disease Control recorded over 8,000 heat-related deaths in the U.S.[^CDCpicture] Heat-related hospitalizations or emergency department visits are estimated to be at least 10 times the death rate.[^Samuelson2020] Heat exhaustion and heat stroke are the most serious heat-related illnesses. Exposure to excessive heat can directly or indirectly cause some illnesses and can exacerbate many preexisting conditions, such as heart and lung disease. People at greatest risk for heat-related illness include children under 5 years of age, people 65 years of age and older, people who are overweight or have existing medical conditions, such as diabetes and heart disease, people who are socially isolated, and low income individuals.[^CDCpicture]
Risk of exposure to excessive heat tends to be higher for people who work outside (e.g., agriculture, construction, and landscaping), and for those living in densely developed urban areas where there is a dearth of vegetation and an abundance of dense materials such as asphalt and concrete that absorb heat and release it more slowly (i.e., urban heat island effect). However, risk of heat-related illness or death varies within urban environments due to a variety of mitigating factors, especially age, race, and wealth.[^Williams2020] Exposure to higher temperatures than a population is accustomed to can also make people more vulnerable to heat-related illnesses and death. Although Massachusetts experiences a generally cool climate, research has shown that people in this state are actually more sensitive to elevated temperatures than those living in warmer states of the country.[^Hondula2015]
Climate change is increasing average global temperatures, as well as the frequency, duration, and severity of extreme heat events.[^IPCCsummary] In the contiguous U.S., annual average temperatures have increased by 1.2°F (0.7°C) over the last few decades and by 1.8°F (1.0°C) relative to the beginning of the 20th century. These changes are happening more quickly in the Northeast. Average annual temperatures have increased by about 3°F (1.7°C) or more in Massachusetts since 1901. By 2050, average annual temperatures in the Northeast are projected to increase by up to 5.1°F (2.8°C) relative to the period 1975–2005, with several more days of extreme heat occurring throughout the state each year.[^4thNatlAssessment] The combination of rapidly increasing temperatures and higher sensitivity puts people in Massachusetts at especially high risk of heat-related illness and death.
The heat risk analysis presented here is based on Land Surface Temperatures (LST) derived from NASA's Moderate-resolution Imaging Spectroradiometer (MODIS) satellite sensor.[^MODIS11] Unlike ambient air temperature measured by ground-based weather stations, LST is a measure of the radiant energy (i.e., emissivity) of the ground or surface. LST values tend to be slightly higher than ambient air temperatures, but are highly correlated with air temperatures. Abundant research has shown that satellite-derived LST is an acceptable proxy for air temperature.[^Kloog2014] Moreover, while weather stations are sparsely and unevenly distributed, satellite-derived LST has the advantage of providing an unbroken and continuous snapshot of temperature at high resolution across the state at any given moment.
The MODIS data used here covers average day and night LST for an 8-day period from July 28 to August 4, 2019 for the entire state. That week constituted the conclusion of a historically warm July for Massachusetts This data originates at 1km spatial resolution and was extracted to average temperatures at each Census Block Group separately for day-night averages, daytime, nighttime, and urban heat island effects. LST data was validated against ground weather station data for the same time periods (see Appendix A).
```{r dataLST, include=FALSE}
library(tmap)
library(tidyverse)
library(sf)
library(tmaptools)
library(maptools)
library(ggcorrplot)
library(spdep)
library(rgdal)
library(gdalUtils)
library(raster)
library(tigris)
options(tigris_use_cache = TRUE, tigris_class = "sf")
library(lwgeom)
library(OpenStreetMap)
# load basic layers
load("DATA/ne_layers.rds")
# # create OSM basemap. see help on openmap for type options.
# ma_bmap <- ne_states_sf_cb %>%
# st_transform(., crs = 4326) %>%
# st_bbox(.) %>%
# read_osm(., type = "stamen-watercolor")
#
# # map basemap
# tm_shape(ma_bmap) + tm_rgb() + tm_shape(ne_states_sf_cb) + tm_borders()
# # STEPS FOR PROCESSING LST DATA
# # MODIS data acquired through USGS EarthExplorer portal. https://earthexplorer.usgs.gov/
# # read in HDF (Hierarchical Data Format) scientific data sets (SDSs) for MODIS from Terra satellite
# # sdsJuly282019h12Terra <- get_subdatasets("DATA/LST/MODISJulyAug2019/MOD11A2.A2019209.h12v04.006.2019218045630.hdf")
# # sdsJuly282019h13Terra <- get_subdatasets("DATA/LST/MODISJulyAug2019/MOD11A2.A2019209.h13v04.006.2019218045706.hdf")
#
# # sdsJuly282019h12 <- get_subdatasets("LST/MODIS8dayLST1km/2345441438/MOD11A2_A2019209_h12v04_006_2019218045630_HEGOUT.hdf")
# # sdsJuly282019h13 <- get_subdatasets("LST/MODIS8dayLST1km/2345441418/MOD11A2_A2019209_h13v04_006_2019218045706_HEGOUT.hdf")
#
# # read in HDF (Hierarchical Data Format) scientific data sets (SDSs) for MODIS from Aqua satellite
# sdsJuly282019h12Aqua <- get_subdatasets("DATA/LST/MODISJulyAug2019/MYD11A2.A2019209.h12v04.006.2019218030935.hdf")
# sdsJuly282019h13Aqua <- get_subdatasets("DATA/LST/MODISJulyAug2019/MYD11A2.A2019209.h13v04.006.2019218033834.hdf")
#
# # retrieve desired data from MODIS Terra (late morning 10:06am - 12:18pm and early evening 8:42pm - 11pm)
# # LSTdayJuly282019h12Terra <- readGDAL(sdsJuly282019h12Terra[1])
# # daytimeH12Terra <- readGDAL(sdsJuly282019h12Terra[3])
# # LSTnightJuly282019h12Terra <- readGDAL(sdsJuly282019h12Terra[5])
# # nighttimeH12Terra <- readGDAL(sdsJuly282019h12Terra[7])
# # LSTdayJuly282019h13Terra <- readGDAL(sdsJuly282019h13Terra[1])
# # daytimeH13Terra <- readGDAL(sdsJuly282019h13Terra[3])
# # LSTnightJuly282019h13Terra <- readGDAL(sdsJuly282019h13Terra[5])
# # nighttimeH13Terra <- readGDAL(sdsJuly282019h13Terra[7])
#
# # retrieve desired data sets from MODIS Aqua (early afternoon 11:48am - 2pm and later evening 12am - 3:06am)
# LSTdayJuly282019h12Aqua <- readGDAL(sdsJuly282019h12Aqua[1])
# daytimeH12Aqua <- readGDAL(sdsJuly282019h12Aqua[3])
# LSTnightJuly282019h12Aqua <- readGDAL(sdsJuly282019h12Aqua[5])
# nighttimeH12Aqua <- readGDAL(sdsJuly282019h12Aqua[7])
# LSTdayJuly282019h13Aqua <- readGDAL(sdsJuly282019h13Aqua[1])
# daytimeH13Aqua <- readGDAL(sdsJuly282019h13Aqua[3])
# LSTnightJuly282019h13Aqua <- readGDAL(sdsJuly282019h13Aqua[5])
# nighttimeH13Aqua <- readGDAL(sdsJuly282019h13Aqua[7])
#
# # assign to raster image to manipulate data for visualization and analysis and convert from Kelvin to Fahrenheit
# LSTdayJuly282019h12r <- raster(LSTdayJuly282019h12Aqua)*(9/5) - 459.67
# LSTnightJuly282019h12r <- raster(LSTnightJuly282019h12Aqua)*(9/5) - 459.67
# LSTdayJuly282019h13r <- raster(LSTdayJuly282019h13Aqua)*(9/5) - 459.67
# LSTnightJuly282019h13r <- raster(LSTnightJuly282019h13Aqua)*(9/5) - 459.67
#
# # # If origins are different, need to adjust tolerance to mosaic
# # origin(LSTdayJuly282019h12r)
# # origin(LSTdayJuly282019h13r)
#
# # mosic rasters
# # create LST day mosaic raster
# LSTdayJuly2019r <- mosaic(LSTdayJuly282019h12r,LSTdayJuly282019h13r,
# fun=mean, tolerance = 1)
#
# # create LST night mosaic raster
# LSTnightJuly2019r <- mosaic(LSTnightJuly282019h12r,LSTnightJuly282019h13r,
# fun=mean, tolerance = 1)
#
# # create average of day and night
# LSTavgJuly2019r <- (LSTdayJuly2019r + LSTnightJuly2019r)/2
#
# # create average day - night difference
# LSTDayNighDiffJuly2019r <- LSTdayJuly2019r - LSTnightJuly2019r
#
# # clean up
# rm(list = ls(pattern = "h12|h13"))
#
# # Crop LST to Massachusetts
# # isolate data for Massachusetts
# ma_blkgrp_sf <- ne_blkgrp_sf %>%
# filter(STATE == "Massachusetts") %>%
# st_transform(., crs = 2805)
#
# ma_state_sf_cb <- ne_states_sf_cb %>%
# filter(NAME == "Massachusetts") %>%
# st_transform(., crs = 2805)
#
# # define new projection
# newproj <- st_crs(ma_blkgrp_sf)[[2]]
#
# LSTavgJuly2019r_ma <- LSTavgJuly2019r %>%
# projectRaster(., crs = newproj) %>%
# crop(., ma_state_sf_cb) %>%
# mask(., ma_state_sf_cb)
#
# LSTdayJuly2019r_ma <- LSTdayJuly2019r %>%
# projectRaster(., crs = newproj) %>%
# crop(., ma_state_sf_cb) %>%
# mask(., ma_state_sf_cb)
#
# LSTnightJuly2019r_ma <- LSTnightJuly2019r %>%
# projectRaster(., crs = newproj) %>%
# crop(., ma_state_sf_cb) %>%
# mask(., ma_state_sf_cb)
#
# LSTDayNighDiffJuly2019r_ma <- LSTDayNighDiffJuly2019r %>%
# projectRaster(., crs = newproj) %>%
# crop(., ma_state_sf_cb) %>%
# mask(., ma_state_sf_cb)
#
# # download Census-designated urban areas to contrast with rural areas
# us_urban_sf <- urban_areas(cb = TRUE)
#
# # create urban area polygon for Massachusetts
# ma_urban_sf <- us_urban_sf %>%
# st_centroid(.) %>%
# st_transform(., crs = 2805) %>%
# st_join(., ma_state_sf_cb, join = st_within, left = FALSE) %>%
# as.data.frame() %>%
# dplyr::select(GEOID10) %>%
# inner_join(us_urban_sf,.,by="GEOID10") %>%
# st_transform(., crs = 2805) %>%
# group_by() %>%
# summarize()
#
# # create rural area polygon
# ma_rural_sf <- ma_state_sf_cb %>%
# group_by() %>%
# summarize() %>%
# st_difference(., ma_urban_sf) %>%
# st_make_valid()
#
# # # select blkgrps that do not intersect wtih urban polygons - rural blkgrps
# # ma_blkgrpLST_rural <- st_disjoint(ma_blkgrp_sf,
# # ma_urban_sf,sparse = FALSE) %>%
# # apply(., 1, all) %>% # convert dense matrix to logical vector
# # ma_blkgrpLST_sf[.,] # subset blkgrps where there is NO intersection
# #
# # # select blkgrps that intersect with urban polygons - urban blkgrps
# # ma_blkgrpLST_urban <- st_intersects(ma_blkgrp_sf,
# # ma_urban_sf,sparse = FALSE) %>%
# # apply(., 1, any) %>% # convert dense matrix to logical vector
# # ma_blkgrpLST_sf[.,] # subset blkgrps where there is ANY intersection
# #
# # ma_blkgrp_urban_simple <- ma_blkgrpLST_urban %>%
# # group_by() %>%
# # summarize()
# #
# # ma_blkgrp_rural_simple <- ma_blkgrpLST_rural %>%
# # group_by() %>%
# # summarize()
#
# # crop LST to urban
# # LSTavgJuly2019r_maUrban <- LSTavgJuly2019r %>%
# # projectRaster(., crs = newproj) %>%
# # crop(., ma_urban_sf) %>%
# # mask(., ma_urban_sf)
# #
# # LSTdayJuly2019r_maUrban <- LSTdayJuly2019r %>%
# # projectRaster(., crs = newproj) %>%
# # crop(., ma_urban_sf) %>%
# # mask(., ma_urban_sf)
# #
# # LSTnightJuly2019r_maUrban <- LSTnightJuly2019r %>%
# # projectRaster(., crs = newproj) %>%
# # crop(., ma_urban_sf) %>%
# # mask(., ma_urban_sf)
# #
# # LSTDayNighDiffJuly2019r_maUrban <- LSTDayNighDiffJuly2019r %>%
# # projectRaster(., crs = newproj) %>%
# # crop(., ma_urban_sf) %>%
# # mask(., ma_urban_sf)
#
# # crop LST to rural
# LSTavgJuly2019r_maRural <- LSTavgJuly2019r %>%
# projectRaster(., crs = newproj) %>%
# crop(., ma_rural_sf) %>%
# mask(., ma_rural_sf)
#
# LSTdayJuly2019r_maRural <- LSTdayJuly2019r %>%
# projectRaster(., crs = newproj) %>%
# crop(., ma_rural_sf) %>%
# mask(., ma_rural_sf)
#
# LSTnightJuly2019r_maRural <- LSTnightJuly2019r %>%
# projectRaster(., crs = newproj) %>%
# crop(., ma_rural_sf) %>%
# mask(., ma_rural_sf)
#
# LSTDayNighDiffJuly2019r_maRural <- LSTDayNighDiffJuly2019r %>%
# projectRaster(., crs = newproj) %>%
# crop(., ma_rural_sf) %>%
# mask(., ma_rural_sf)
#
# # Identify stats to compute thresholds
# # LSTavgJuly2019mean <- cellStats(x = LSTavgJuly2019r_ne, stat = "mean")
# # LSTavgJuly2019Urbanmean <- cellStats(x = LSTavgJuly2019r_neUrban, stat = "mean")
# LSTavgJuly2019Ruralmean <- cellStats(x = LSTavgJuly2019r_maRural, stat = "mean")
# # LSTdayJuly2019mean <- cellStats(x = LSTdayJuly2019r_ne, stat = "mean")
# # LSTdayJuly2019Urbanmean <- cellStats(x = LSTdayJuly2019r_neUrban, stat = "mean")
# LSTdayJuly2019Ruralmean <- cellStats(x = LSTdayJuly2019r_maRural, stat = "mean")
# # LSTnightJuly2019mean <- cellStats(x = LSTnightJuly2019r_ne, stat = "mean")
# # LSTnightJuly2019Urbanmean <- cellStats(x = LSTnightJuly2019r_neUrban, stat = "mean")
# LSTnightJuly2019Ruralmean <- cellStats(x = LSTnightJuly2019r_maRural, stat = "mean")
# # LSTavgJuly2019sd <- cellStats(x = LSTavgJuly2019r_ne, stat = "sd")
# # LSTavgJuly2019Urbansd <- cellStats(x = LSTavgJuly2019r_neUrban, stat = "sd")
# # LSTavgJuly2019Ruralsd <- cellStats(x = LSTavgJuly2019r_neRural, stat = "sd")
# # LSTdayJuly2019sd <- cellStats(x = LSTdayJuly2019r_ne, stat = "sd")
# # LSTdayJuly2019Urbansd <- cellStats(x = LSTdayJuly2019r_neUrban, stat = "sd")
# # LSTdayJuly2019Ruralsd <- cellStats(x = LSTdayJuly2019r_neRural, stat = "sd")
# # LSTnightJuly2019sd <- cellStats(x = LSTnightJuly2019r_ne, stat = "sd")
# # LSTnightJuly2019Urbansd <- cellStats(x = LSTnightJuly2019r_neUrban, stat = "sd")
# # LSTnightJuly2019Ruralsd <- cellStats(x = LSTnightJuly2019r_neRural, stat = "sd")
#
# # # Create rasters of cells exceeding threshold
# # LSTavgJulyHI <- LSTavgJuly2019r_ne >
# # (LSTavgJuly2019mean + LSTavgJuly2019sd * 3)
# # LSTdayJulyHI <- LSTdayJuly2019r_ne >
# # (LSTdayJuly2019mean + LSTdayJuly2019sd * 3)
# # LSTnightJulyHI <- LSTnightJuly2019r_ne >
# # (LSTnightJuly2019mean + LSTnightJuly2019sd * 3)
#
# # Extract raster values to block groups
# # check for empty geometries
# any(is.na(st_dimension(ma_blkgrp_sf)))
# # identiy empty geometries
# empty_geo <- st_is_empty(ma_blkgrp_sf)
# # filter out empty geometries
# ma_blkgrp_sf <- ma_blkgrp_sf[!empty_geo,]
# # clean up
# rm(empty_geo)
#
# # Extract mean LST values within each block group to a df
# meanLSTavg_df <- ma_blkgrp_sf %>%
# dplyr::select(GEOID) %>%
# as_Spatial() %>%
# extract(LSTavgJuly2019r_ma, .,
# fun=mean, sp=TRUE, na.rm=TRUE, small=TRUE) %>%
# as.data.frame() %>%
# rename(meanAvgLST = layer)
#
# meanLSTday_df <- ma_blkgrp_sf %>%
# dplyr::select(GEOID) %>%
# as_Spatial() %>%
# extract(LSTdayJuly2019r_ma, .,
# fun=mean, sp=TRUE, na.rm=TRUE, small=TRUE) %>%
# as.data.frame() %>%
# rename(meanDayLST = layer)
#
# meanLSTnight_df <- ma_blkgrp_sf %>%
# dplyr::select(GEOID) %>%
# as_Spatial() %>%
# extract(LSTnightJuly2019r_ma, .,
# fun=mean, sp=TRUE, na.rm=TRUE, small=TRUE) %>%
# as.data.frame() %>%
# rename(meanNightLST = layer)
#
# meanLSTDayNightDiff_df <- ma_blkgrp_sf %>%
# dplyr::select(GEOID) %>%
# as_Spatial() %>%
# extract(LSTDayNighDiffJuly2019r_ma,.,
# fun=mean, sp=TRUE, na.rm=TRUE, small=TRUE) %>%
# as.data.frame() %>%
# rename(meanDayNightDiffLST = layer)
#
# # # Extract mean LST value for rural areas in order to compute urban heat island
# # meanLSTavgRural_df <- ma_blkgrp_rural_simple %>%
# # as_Spatial() %>%
# # extract(LSTavgJuly2019r_ne, .,
# # fun=mean, sp=TRUE, na.rm=TRUE, small=TRUE) %>%
# # as.data.frame() %>%
# # rename(meanAvgLSTrural = layer)
#
# # Join LST statistics to block groups
# ma_blkgrpLST_sf <- ma_blkgrp_sf %>%
# left_join(., meanLSTavg_df, by = "GEOID") %>%
# left_join(., meanLSTday_df, by = "GEOID") %>%
# left_join(., meanLSTnight_df, by = "GEOID") %>%
# left_join(., meanLSTDayNightDiff_df, by = "GEOID") %>%
# mutate(UHI24avg = meanAvgLST - LSTavgJuly2019Ruralmean,
# UHIDay = meanDayLST - LSTdayJuly2019Ruralmean,
# UHINight = meanNightLST - LSTnightJuly2019Ruralmean)
#
# save(ma_blkgrpLST_sf, ma_urban_sf, ma_rural_sf, file = "DATA/LST/ma_blkgrpLST_sf.rds")
#
# # clean up rasters
# rm(list = ls(pattern = "July2019r"))
#
# load previously processed LST data
load("DATA/LST/ma_blkgrpLST_sf.rds")
ma_towns_sf_pts <- county_subdivisions(state = "MA", cb = TRUE) %>%
filter(NAME %in% c("Boston",
"Salem",
"Lawrence",
"Lowell",
"Newburyport",
"Rockport",
"Brockton",
"New Bedford",
"Plymouth",
"Worcester",
"Springfield",
"Pittsfield",
"Athol",
"Wrentham,",
"Great Barrington",
"Holyoke",
"Barnstable Town")) %>%
st_transform(., crs = 2805) %>%
st_centroid(of_largest_polygon = TRUE)
# Create road layer for context
ma_highways <- primary_roads() %>%
filter(FULLNAME %in% c("I- 84","I- 90","I- 91","I- 95","I- 190","I- 195","I- 290","I- 395","I- 495","US Hwy 6","US Hwy 202","Mohawk Trl","George W Stanton Hwy","State Rte 2","Mass State Hwy","Concord Tpke","State Rte 25")) %>%
tmaptools::crop_shape(., ne_states_sf_cb) %>%
st_transform(., crs = 2805)
ma_highways2nd <- primary_secondary_roads("MA") %>%
filter(FULLNAME %in% c("US Hwy 6","Mohawk Trl","State Rte 2","Cambridge Tpke")) %>%
st_transform(., crs = 2805)
# Extract highway segments for labeling
I90roadSegment <- ma_highways %>%
filter(LINEARID == "1103745154991")
I90roadSegment2 <- ma_highways %>%
filter(LINEARID == "110340769311")
I91roadSegment <- ma_highways %>%
filter(LINEARID == "1104748241453")
I95roadSegment <- ma_highways %>%
filter(LINEARID == "1105569136116")
I95roadSegment2 <- ma_highways %>%
filter(LINEARID == "1103737956638")
I195roadSegment2 <- ma_highways %>%
filter(LINEARID == "1101922014382")
I395roadSegment <- ma_highways %>%
filter(LINEARID == "1104259933162")
I495roadSegment <- ma_highways %>%
filter(LINEARID == "1103745404033")
I495roadSegment2 <- ma_highways %>%
filter(LINEARID == "1105589457557")
I495roadSegment3 <- ma_highways %>%
filter(LINEARID == "1101922014436")
StRt2Segment <- ma_highways2nd %>%
filter(LINEARID == "1106087431756")
USHwy6Segment <- ma_highways %>%
filter(LINEARID == "1109096413415")
# Create custom icons of highway shields
I90 <- tmap_icons(file = "https://upload.wikimedia.org/wikipedia/commons/thumb/c/ca/I-90.svg/200px-I-90.svg.png")
I95 <- tmap_icons(file = "https://upload.wikimedia.org/wikipedia/commons/thumb/6/61/I-95.svg/200px-I-95.svg.png")
I395 <- tmap_icons("https://upload.wikimedia.org/wikipedia/commons/thumb/1/16/I-395.svg/200px-I-395.svg.png")
I91 <- tmap_icons("https://upload.wikimedia.org/wikipedia/commons/thumb/9/90/I-91.svg/200px-I-91.svg.png")
I495 <- tmap_icons("https://upload.wikimedia.org/wikipedia/commons/thumb/8/8c/I-495.svg/200px-I-495.svg.png")
Hwy2 <- tmap_icons("https://upload.wikimedia.org/wikipedia/commons/thumb/5/54/MA_Route_2.svg/240px-MA_Route_2.svg.png")
Hwy6 <- tmap_icons("https://upload.wikimedia.org/wikipedia/commons/thumb/e/ef/US_6.svg/200px-US_6.svg.png")
# Calculate percentiles for pops of concern and temperature and create ranks for comparison
ma_blkgrpLST_sf <- ma_blkgrpLST_sf %>%
mutate(pctile24hrLST = percent_rank(meanAvgLST),
rank24hrLST = case_when(
pctile24hrLST < 0.2 ~ 1,
pctile24hrLST >= 0.2 & pctile24hrLST < 0.4 ~ 2,
pctile24hrLST >= 0.4 & pctile24hrLST < 0.6 ~ 3,
pctile24hrLST >= 0.6 & pctile24hrLST < 0.8 ~ 4,
pctile24hrLST >= 0.8 ~ 5),
pctileDayLST = percent_rank(meanDayLST),
rankDayLST = case_when(
pctileDayLST < 0.2 ~ 1,
pctileDayLST >= 0.2 & pctileDayLST < 0.4 ~ 2,
pctileDayLST >= 0.4 & pctileDayLST < 0.6 ~ 3,
pctileDayLST >= 0.6 & pctileDayLST < 0.8 ~ 4,
pctileDayLST >= 0.8 ~ 5),
pctileNightLST = percent_rank(meanNightLST),
rankNightLST = case_when(
pctileNightLST < 0.2 ~ 1,
pctileNightLST >= 0.2 & pctileNightLST < 0.4 ~ 2,
pctileNightLST >= 0.4 & pctileNightLST < 0.6 ~ 3,
pctileNightLST >= 0.6 & pctileNightLST < 0.8 ~ 4,
pctileNightLST >= 0.8 ~ 5),
pctileMinority = percent_rank(minority_pctE),
rankMinority = case_when(
pctileMinority < 0.2 ~ 1,
pctileMinority >= 0.2 & pctileMinority < 0.4 ~ 2,
pctileMinority >= 0.4 & pctileMinority < 0.6 ~ 3,
pctileMinority >= 0.6 & pctileMinority < 0.8 ~ 4,
pctileMinority >= 0.8 ~ 5),
LST24MinorScore = (rank24hrLST + rankMinority)/2,
LSTDayMinorScore = (rankDayLST + rankMinority)/2,
LSTNightMinorScore = (rankNightLST + rankMinority)/2,
pctileEngLang = percent_rank(eng_limit_pctE),
rankEngLang = case_when(
pctileEngLang < 0.2 ~ 1,
pctileEngLang >= 0.2 & pctileEngLang < 0.4 ~ 2,
pctileEngLang >= 0.4 & pctileEngLang < 0.6 ~ 3,
pctileEngLang >= 0.6 & pctileEngLang < 0.8 ~ 4,
pctileEngLang >= 0.8 ~ 5),
LST24EngLangScore = (rank24hrLST + rankEngLang)/2,
LSTDayEngLangScore = (rankDayLST + rankEngLang)/2,
LSTNightEngLangScore = (rankNightLST + rankEngLang)/2,
pctileLowInc = percent_rank(pct2povE),
rankLowInc = case_when(
pctileLowInc < 0.2 ~ 1,
pctileLowInc >= 0.2 & pctileLowInc < 0.4 ~ 2,
pctileLowInc >= 0.4 & pctileLowInc < 0.6 ~ 3,
pctileLowInc >= 0.6 & pctileLowInc < 0.8 ~ 4,
pctileLowInc >= 0.8 ~ 5),
LST24LowIncScore = (rank24hrLST + rankLowInc)/2,
LSTDayLowIncScore = (rankDayLST + rankLowInc)/2,
LSTNightLowIncScore = (rankNightLST + rankLowInc)/2,
pctileNoHSDip = percent_rank(pct_lthsE),
rankNoHSDip = case_when(
pctileNoHSDip < 0.2 ~ 1,
pctileNoHSDip >= 0.2 & pctileNoHSDip < 0.4 ~ 2,
pctileNoHSDip >= 0.4 & pctileNoHSDip < 0.6 ~ 3,
pctileNoHSDip >= 0.6 & pctileNoHSDip < 0.8 ~ 4,
pctileNoHSDip >= 0.8 ~ 5),
LST24NoHSDipScore = (rank24hrLST + rankNoHSDip)/2,
LSTDayNoHSDipScore = (rankDayLST + rankNoHSDip)/2,
LSTNightNoHSDipScore = (rankNightLST + rankNoHSDip)/2,
pctileUnder5 = percent_rank(pct_under5E),
rankUnder5 = case_when(
pctileUnder5 < 0.2 ~ 1,
pctileUnder5 >= 0.2 & pctileUnder5 < 0.4 ~ 2,
pctileUnder5 >= 0.4 & pctileUnder5 < 0.6 ~ 3,
pctileUnder5 >= 0.6 & pctileUnder5 < 0.8 ~ 4,
pctileUnder5 >= 0.8 ~ 5),
LST24Under5Score = (rank24hrLST + rankUnder5)/2,
LSTDayUnder5Score = (rankDayLST + rankUnder5)/2,
LSTNightUnder5Score = (rankNightLST + rankUnder5)/2,
pctileOver64 = percent_rank(pct_over64E),
rankOver64 = case_when(
pctileOver64 < 0.2 ~ 1,
pctileOver64 >= 0.2 & pctileOver64 < 0.4 ~ 2,
pctileOver64 >= 0.4 & pctileOver64 < 0.6 ~ 3,
pctileOver64 >= 0.6 & pctileOver64 < 0.8 ~ 4,
pctileOver64 >= 0.8 ~ 5),
LST24Over64Score = (rank24hrLST + rankOver64)/2,
LSTDayOver64Score = (rankDayLST + rankOver64)/2,
LSTNightOver64Score = (rankNightLST + rankOver64)/2
)
# Calculate totals of priority populations for use in computing percentages in later tables
statepopsdf <- ma_blkgrpLST_sf %>%
as.data.frame() %>%
mutate(MA_LOWINC = if_else(MA_INCOME == "I", totalpopE, 0)) %>%
mutate(MA_LOWINC = replace_na(MA_LOWINC,0)) %>%
mutate(MA_MINORITY = if_else(MA_MINORITY == "M", totalpopE,0)) %>%
mutate(MA_MINORITY = replace_na(MA_MINORITY,0)) %>%
mutate(MA_ENGLISH = if_else(MA_ENGLISH == "E", totalpopE,0)) %>%
mutate(MA_ENGLISH = replace_na(MA_ENGLISH,0)) %>%
group_by(STATE) %>%
summarize(TotalPop = sum(totalpopE, na.rm = TRUE),
TotalHH = sum(householdsE, na.rm = TRUE),
TotalLEH = sum(eng_limitE, na.rm = TRUE),
TotalMin = sum(minorityE, na.rm = TRUE),
TotalLowInc = sum(num2povE, na.rm = TRUE),
TotalNoHS = sum(lthsE, na.rm = TRUE),
TotalOver64 = sum(over64E, na.rm = TRUE),
TotalUnder5 = sum(under5E, na.rm = TRUE),
TotalMA_LOWINC = sum(MA_LOWINC, na.rm = TRUE),
TotalMA_MINORITY = sum(MA_MINORITY, na.rm = TRUE),
TotalMA_ENGLISH = sum(MA_ENGLISH, na.rm = TRUE))
# ma_towns_sf <- ne_towns_sf %>%
# dplyr::select(GEOID,NAME) %>%
# st_transform(., crs = 2805)
ma_towns_sf <- county_subdivisions(state = "MA", cb = TRUE) %>%
st_transform(., crs = 2805)
townpopsdf <- ma_blkgrpLST_sf %>%
dplyr::select(-GEOID) %>%
st_centroid(.) %>%
st_intersection(ma_towns_sf,.) %>%
as.data.frame() %>%
mutate(MA_LOWINC = if_else(MA_INCOME == "I", totalpopE, 0)) %>%
mutate(MA_LOWINC = replace_na(MA_LOWINC,0)) %>%
mutate(MA_MINORITY = if_else(MA_MINORITY == "M", totalpopE,0)) %>%
mutate(MA_MINORITY = replace_na(MA_MINORITY,0)) %>%
mutate(MA_ENGLISH = if_else(MA_ENGLISH == "E", totalpopE,0)) %>%
mutate(MA_ENGLISH = replace_na(MA_ENGLISH,0)) %>%
group_by(NAME) %>%
summarize(GEOID = unique(GEOID),
TotalPop = sum(totalpopE, na.rm = TRUE),
TotalHH = sum(householdsE, na.rm = TRUE),
TotalLEH = sum(eng_limitE, na.rm = TRUE),
TotalMin = sum(minorityE, na.rm = TRUE),
TotalLowInc = sum(num2povE, na.rm = TRUE),
TotalNoHS = sum(lthsE, na.rm = TRUE),
TotalOver64 = sum(over64E, na.rm = TRUE),
TotalUnder5 = sum(under5E, na.rm = TRUE),
TotalMA_LOWINC = sum(MA_LOWINC, na.rm = TRUE),
TotalMA_MINORITY = sum(MA_MINORITY, na.rm = TRUE),
TotalMA_ENGLISH = sum(MA_ENGLISH, na.rm = TRUE))
ma_state_sf_cb <- ne_states_sf_cb %>%
filter(NAME == "Massachusetts") %>%
st_transform(., crs = st_crs(ma_blkgrpLST_sf))
```
## Day-Night Average Land Surface Temperatures
Day-Night Average Land Surface Temperatures (LST) represent a simple average of LST values collected during the day (11:48am - 2pm) and during the night (12am - 3:06am) over eight days, from July 28 to August 4, 2019. These day-night average temperatures varied significantly across Massachusetts, with highest temperatures apparent in the most urbanized areas of the state. This is particularly apparent in the greater Boston area, and also around the state's largest cities, including Springfield, Worcester, Lowell, Lawrence, Brockton, and New Bedford (see Figure \@ref(fig:mapLSTavgMA)).
```{r mapLSTavgMA, fig.align = "center", fig.cap="Map of July-Aug 2019 day-night Average Land Surface Temperatures (LST) across Massachusetts at Census Block Group level.", strip.white=TRUE}
# Map of average LST across Massachusetts
m <- tm_shape(ma_blkgrpLST_sf, unit = "mi") +
tm_fill("meanAvgLST", style = "quantile", palette = "Oranges",
title = expression("Temperature " ( degree*F)),
legend.hist = FALSE,
colorNA = NULL,
textNA = NULL,
legend.format=list(list(digits=1)),
legend.is.portrait = TRUE) +
tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.2, alpha = 0.8) +
tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
tm_shape(ma_highways) + tm_lines(col = "seashell4", lwd = 1) +
tm_shape(ma_highways2nd) + tm_lines(col = "seashell4", lwd = 1) +
tm_shape(I95roadSegment) +
tm_symbols(shape = I95, border.lwd = NA, size = .1) +
tm_shape(I95roadSegment2) +
tm_symbols(shape = I95, border.lwd = NA, size = .1) +
tm_shape(I395roadSegment) +
tm_symbols(shape = I395, border.lwd = NA, size = .1) +
tm_shape(I91roadSegment) +
tm_symbols(shape = I91, border.lwd = NA, size = .1) +
tm_shape(I495roadSegment) +
tm_symbols(shape = I495, border.lwd = NA, size = .1) +
tm_shape(I495roadSegment2) +
tm_symbols(shape = I495, border.lwd = NA, size = .1) +
tm_shape(I495roadSegment3) +
tm_symbols(shape = I495, border.lwd = NA, size = .1) +
tm_shape(I90roadSegment) +
tm_symbols(shape = I90, border.lwd = NA, size = .1) +
tm_shape(I90roadSegment2) +
tm_symbols(shape = I90, border.lwd = NA, size = .1) +
tm_shape(ma_towns_sf_pts) + tm_dots() +
tm_text("NAME", size = 0.5, col = "black",
xmod = 0.7, ymod = 0.2, shadow = TRUE) +
tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
position = c(0.6,0.005)) +
tm_layout(title = "Day-Night Avg\nLand Surface\nTemperatures\nJuly-Aug 2019",
frame = FALSE, main.title.size = 0.8,
legend.outside = TRUE,
legend.title.size = 0.8,
legend.outside.position = c("right", "top"),
legend.hist.width = 0.9)
tmap_save(m, "DATA/LST/ma_24hrLST.png",
height = 4, width = 8, units = "in", dpi = 600)
knitr::include_graphics("DATA/LST/ma_24hrLST.png")
```
```{r LST24outliers, include=FALSE}
# Identify municipalities with high heat outliers in Massachusetts
ma_LST24outliers <- ma_blkgrpLST_sf %>%
arrange(desc(meanAvgLST)) %>%
slice(1:10) %>%
st_centroid(.) %>%
st_intersection(ma_towns_sf,.) %>%
as.data.frame() %>%
group_by(NAME) %>%
summarize(unique(NAME)) %>%
pull(NAME) %>%
paste(.,sep = " ", collapse = ", ")
```
Figure \@ref(fig:boxplotLSTavgMA) is a boxplot of average day-night temperatures by Block Group . The box represents temperatures ranging between the 25th and 75th percentiles. The line that divides the box into 2 parts represents the median temperature for all Block Groups, which in this case is `r round(median(ma_blkgrpLST_sf$meanAvgLST, na.rm = TRUE),2)`°. Half of the state's Block Groups are below the median and half are above the median. Each dot represents an individual Block Group. The top ten highest outlier Block Groups (far right) were found in the following municipalities: `r ma_LST24outliers`.
```{r boxplotLSTavgMA, fig.align = "center", fig.cap="Boxplot of July-Aug 2019 day-night average Land Surface Temperatures (LST) across Massachusetts at Census Block Group level."}
# boxplot of PM2.5 by state for 2016
ma_blkgrpLST_sf %>%
as.data.frame() %>%
drop_na(meanAvgLST) %>%
ggplot(aes(x = STATE, y = meanAvgLST, fill = STATE)) +
geom_boxplot() +
geom_jitter(color = "black", size = 0.4, alpha = 0.6) +
ggtitle("Day-Night Average Land Surface Temperatures July-Aug 2019") +
theme_minimal() +
scale_y_continuous(labels = function(x) paste0(x, "°")) +
theme(legend.position = "none", axis.text=element_text(size=8)) +
xlab(NULL) +
ylab("Temperature (°F)") +
coord_flip()
```
### Day-Night Average Urban Heat Island Effect
The urban heat island (UHI) effect is a phenomenon in which temperatures in urban areas tend to be higher than surrounding non-urban or rural areas. These elevated urban air temperatures are a consequence of the high density of buildings, roads, and other impervious infrastructure in urban areas that absorb heat and release it more slowly. The UHI effect is compounded by the relative absence of vegetation which would otherwise moderate temperatures through evaporative cooling.[^Zhou2019] Higher summer temperatures due to UHI increase the risk of heat-related illnesses, result in increased energy use for air conditioning, and exacerbate air pollution, especially ground-level ozone. As the climate warms, UHI are likely to exacerbate both temperature and air pollution risks.
The UHI effect analyzed here is defined as the difference in temperature between urbanized areas and rural areas in Massachusetts.[^urban] Specifically, the average temperature of the rural areas of Massachusetts is subtracted from the average temperatures of each Census Block Group. The resulting value shows how much warmer or cooler each Census Block Group is when compared to the rural background reference. Positive values indicate that a Census Block Group is warmer than the rural background average and may be experiencing the UHI effect. The larger the value, the stronger the effect.
The UHI effect for day-night average temperatures in Massachusetts are most pronounced in the inner core of urbanized areas (see Figure \@ref(fig:mapUHIavgMA)), particularly around Boston and extending north to Lowell and Lawrence, and south toward Brockton, Fall River, and New Bedford. Temperatures in these areas are up to `r round(max(ma_blkgrpLST_sf$UHI24avg, na.rm = TRUE),1)`° warmer than the rural background average. By contrast, most of the western half of the state, except for the Springfield area, is closer to the rural average.
```{r mapUHIavgMA, fig.align = "center", fig.cap="Map of July-Aug 2019 day-night average Urban Heat Island (UHI) effect across Massachusetts at Census Block Group level. UHI is the difference between Block Group average temperature and average temperature of rural areas.", strip.white=TRUE}
# Map of average LST UHI across Massachusetts
m <- tm_shape(ma_blkgrpLST_sf, unit = "mi") +
tm_fill("UHI24avg", palette = "-RdYlGn", midpoint = 0,
legend.reverse = TRUE,
title = expression("Temperature " ( degree*F)),
legend.hist = FALSE,
colorNA = NULL,
textNA = NULL,
legend.format=list(list(digits=1)),
legend.is.portrait = TRUE) +
tm_shape(ne_states_sf_cb) + tm_borders(lwd = 0.2, alpha = 0.8) +
tm_text("STUSPS", size = 0.7, remove.overlap = TRUE, col = "gray") +
tm_shape(ma_highways) + tm_lines(col = "seashell4", lwd = 1) +
tm_shape(ma_highways2nd) + tm_lines(col = "seashell4", lwd = 1) +
tm_shape(I95roadSegment) +
tm_symbols(shape = I95, border.lwd = NA, size = .1) +
tm_shape(I95roadSegment2) +
tm_symbols(shape = I95, border.lwd = NA, size = .1) +
tm_shape(I395roadSegment) +
tm_symbols(shape = I395, border.lwd = NA, size = .1) +
tm_shape(I91roadSegment) +
tm_symbols(shape = I91, border.lwd = NA, size = .1) +
tm_shape(I495roadSegment) +
tm_symbols(shape = I495, border.lwd = NA, size = .1) +
tm_shape(I495roadSegment2) +
tm_symbols(shape = I495, border.lwd = NA, size = .1) +
tm_shape(I495roadSegment3) +
tm_symbols(shape = I495, border.lwd = NA, size = .1) +
tm_shape(I90roadSegment) +
tm_symbols(shape = I90, border.lwd = NA, size = .1) +
tm_shape(I90roadSegment2) +
tm_symbols(shape = I90, border.lwd = NA, size = .1) +
tm_shape(ma_towns_sf_pts) + tm_dots() +
tm_text("NAME", size = 0.5, col = "black",
xmod = 0.7, ymod = 0.2, shadow = TRUE) +
tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.5,
position = c(0.6,0.005)) +
tm_layout(title = "Day-Night Avg\nUrban Heat\nIsland\nJuly-Aug 2019",
frame = FALSE, main.title.size = 0.8,
legend.outside = TRUE,
legend.title.size = 0.8,
legend.outside.position = c("right", "top"),
legend.hist.width = 0.9)
tmap_save(m, "DATA/LST/ma_24hrUHI.png",
height = 4, width = 8, units = "in", dpi = 600)
knitr::include_graphics("DATA/LST/ma_24hrUHI.png")
```
```{r LST24UHIoutliers, include=FALSE}
ma_LST24UHIoutliers <- ma_blkgrpLST_sf %>%
arrange(desc(UHI24avg)) %>%
slice(1:10) %>%
st_centroid(.) %>%
st_intersection(ma_towns_sf,.) %>%
as.data.frame() %>%
group_by(NAME) %>%
summarize(unique(NAME)) %>%
pull(NAME) %>%
paste(.,sep = " ", collapse = ", ")
```
Figure \@ref(fig:boxplotUHIavgMA) is a boxplot of day-night average Urban Heat Island (UHI) effect by Block Group. The box represents UHI values ranging between the 25th and 75th percentiles. The line that divides the box into 2 parts represents the median UHI for all Block Groups, which in this case is `r round(median(ma_blkgrpLST_sf$UHI24avg, na.rm = TRUE),2)`°. Half of the state's Block Groups are below the median and half are above the median. Each dot represents an individual Block Group. The top ten highest outlier Block Groups (far right) were found in the following municipalities: `r ma_LST24UHIoutliers`.
```{r boxplotUHIavgMA, fig.align = "center", fig.cap="Boxplot of July-Aug 2019 day-night average Urban Heat Island (UHI) across Massachusetts at Census Block Group level. UHI is the difference between Block Group average temperature and average temperature of rural areas."}
# boxplot of PM2.5 by state for 2016
ma_blkgrpLST_sf %>%
as.data.frame() %>%
drop_na(UHI24avg) %>%
ggplot(aes(x = STATE, y = UHI24avg, fill = STATE)) +
geom_boxplot() +
geom_jitter(color = "black", size = 0.4, alpha = 0.6) +
ggtitle("Day-Night Average UHI July-Aug 2019") +
theme_minimal() +
scale_y_continuous(labels = function(x) paste0(x, "°")) +
theme(legend.position = "none", axis.text=element_text(size=8)) +
xlab(NULL) +
ylab("Temperature (°F)") +
coord_flip()
```
### Day-Night Average Land Surface Temperatures and Priority Populations
```{r LST24hrTempPop, include=FALSE}
# Pop Weighted avg of 24hr LST exposure for all Groups in Massachusetts relative to MA average
LST24hrPop <- ma_blkgrpLST_sf %>%
as.data.frame() %>%
mutate(MA_LOWINC = if_else(MA_INCOME == "I", totalpopE, 0)) %>%
mutate(MA_LOWINC = replace_na(MA_LOWINC,0)) %>%
mutate(MA_MINORITY = if_else(MA_MINORITY == "M", totalpopE,0)) %>%
mutate(MA_MINORITY = replace_na(MA_MINORITY,0)) %>%
mutate(MA_ENGLISH = if_else(MA_ENGLISH == "E", totalpopE,0)) %>%
mutate(MA_ENGLISH = replace_na(MA_ENGLISH,0)) %>%
# filter(STATE == "Massachusetts") %>%
dplyr::select(totalpopE,
nhwhitepopE,
minorityE,
nhblackpopE,
nhamerindpopE,
nhasianpopE,
nhnativhpopE,
nhotherpopE,
nh2morepopE,
hisppopE,
povknownE,
num2povE,
eng_hhE,
eng_limitE,
age25upE,
lthsE,
allAgesE,
under5E,
over64E,
MA_LOWINC,
MA_MINORITY,
MA_ENGLISH,
meanAvgLST) %>%
gather(key = Group, value = Pop, totalpopE:MA_ENGLISH) %>%
group_by(Group) %>%
summarize(LST24hrwMean = weighted.mean(x = meanAvgLST, w = Pop, na.rm = TRUE),
LST24hrMean = mean(meanAvgLST, na.rm = TRUE)) %>%
filter(Group %in% c("minorityE","eng_limitE","num2povE","lthsE","under5E","over64E","MA_LOWINC","MA_MINORITY","MA_ENGLISH")) %>%
mutate(Group = case_when(
Group == "minorityE" ~ "Minority",
Group == "eng_limitE" ~ "Limited English HH",
Group == "num2povE" ~ "Low Income",
Group == "lthsE" ~ "No HS Dip",
Group == "under5E" ~ "Under 5",
Group == "over64E" ~ "Over 64",
Group == "MA_LOWINC" ~ "MA Low Income",
Group == "MA_MINORITY" ~ "MA Minority",
Group == "MA_ENGLISH" ~ "MA Limited English HH"))
# Pop Weighted avg of LST24hr for all Groups in Massachusetts relative to MA average
LST24hrPopTab <- ma_blkgrpLST_sf %>%
as.data.frame() %>%
mutate(MA_LOWINC = if_else(MA_INCOME == "I", totalpopE, 0)) %>%
mutate(MA_LOWINC = replace_na(MA_LOWINC,0)) %>%
mutate(MA_MINORITY = if_else(MA_MINORITY == "M", totalpopE,0)) %>%
mutate(MA_MINORITY = replace_na(MA_MINORITY,0)) %>%
mutate(MA_ENGLISH = if_else(MA_ENGLISH == "E", totalpopE,0)) %>%
mutate(MA_ENGLISH = replace_na(MA_ENGLISH,0)) %>%
# filter(STATE == "Massachusetts") %>%
dplyr::select(totalpopE,
nhwhitepopE,
minorityE,
nhblackpopE,
nhamerindpopE,
nhasianpopE,
nhnativhpopE,
nhotherpopE,
nh2morepopE,
hisppopE,
povknownE,
num2povE,
eng_hhE,
eng_limitE,
age25upE,
lthsE,
allAgesE,
under5E,
over64E,
MA_LOWINC,
MA_MINORITY,
MA_ENGLISH,
meanAvgLST) %>%
gather(key = Group, value = Pop, totalpopE:MA_ENGLISH) %>%
group_by(Group) %>%
summarize(LST24hrwMean = weighted.mean(x = meanAvgLST, w = Pop, na.rm = TRUE),
LST24hrMean = mean(meanAvgLST, na.rm = TRUE)) %>%
spread(key = Group, value = LST24hrwMean) %>%
transmute(Minority = (minorityE/LST24hrMean - 1)*100,
#Minority_NHW = (minorityE/nhwhitepopE - 1)*100,
`Limited English HH` = (eng_limitE/LST24hrMean - 1)*100,
`Low Income` = (num2povE/LST24hrMean - 1)*100,
`No HS Dip` = (lthsE/LST24hrMean - 1)*100,
`Under 5` = (under5E/LST24hrMean - 1)*100,
`Over 64` = (over64E/LST24hrMean - 1)*100,
`MA Low Income` = (MA_LOWINC/LST24hrMean - 1)*100,
`MA Minority` = (MA_MINORITY/LST24hrMean - 1)*100,
`MA Limited English HH` = (MA_ENGLISH/LST24hrMean - 1)*100) %>%
gather(key = Group, value = Pct) %>%
left_join(x = LST24hrPop, y = ., by = "Group") %>%
mutate(Difference = LST24hrwMean - LST24hrMean) %>%
transmute(Group = Group,
LST24hrwMean = paste0(round(LST24hrwMean,2),"°"),
Difference = if_else(Difference > 0,
paste0("+",round(Difference,2),"°"),
paste0(round(Difference,2),"°")),
Pct = if_else(Pct > 0, paste0("+",round(Pct,2),"%"),
paste0(round(Pct,2),"%")))
```
In addition to variations in the general geography of Land Surface Temperatures (LST), temperature exposure also varies demographically. Figure \@ref(fig:plotLST24hrTempPopAvgMA) and Table \@ref(tab:tabLST24hrPopAvgMA) show population-weighted exposures for priority populations relative to the day-night average LST for the state of `r round(mean(ma_blkgrpLST_sf$meanAvgLST,na.rm=TRUE),1)`°. For example, limited English speaking households in Massachusetts, as defined by the state's environmental justice policy, were exposed to a population-weighted day-night average LST of `r round(max(LST24hrPop$LST24hrwMean),1)`°, approximately `r round(max(LST24hrPop$LST24hrwMean) - max(LST24hrPop$LST24hrMean),1)`° (`r LST24hrPopTab[3,4]`) above the state day-night average LST. Similarly, People of Color, as defined by the state's environmental justice policy, were exposed to a population-weighted day-night average LST of `r LST24hrPop %>% filter(Group == "MA Minority") %>% summarize(round(max(LST24hrwMean),1)) %>% pull()`°, over `r round(LST24hrPop %>% filter(Group == "MA Minority") %>% summarize(max(LST24hrwMean)) %>% pull() - max(LST24hrPop$LST24hrMean),1)`° (`r LST24hrPopTab[5,4]`) above the state day-night average LST. By contrast, persons over age 64 were, on average, exposed to population-weighted day-night average LSTs of `r LST24hrPopTab[8,3]` below the state average.
```{r plotLST24hrTempPopAvgMA, echo=FALSE, message=FALSE, warning=FALSE, fig.align = "center", fig.cap="Population-weighted average exposures to day-night average temperatures for priority populations in Massachusetts relative to the state average."}
LST24hrPop %>%
mutate(Group = recode(Group, "Minority" = "People of Color",
"MA Minority" = "MA POC",
"No HS Dip" = "No HS Diploma")) %>%
ggplot(aes(x = reorder(Group,LST24hrwMean),
y = LST24hrwMean)) +
geom_segment(aes(x = reorder(Group,LST24hrwMean),
xend = reorder(Group,LST24hrwMean),
y = LST24hrMean, yend = LST24hrwMean),
color = "firebrick1") +
geom_point(color = "firebrick4", size = 4, alpha = 0.8) +
coord_flip() +
labs(x = "", y = "", title = "Population-Weighted Temperature Exposure to\nDay-Night Average LST (°F)") +
theme_light() +
theme(panel.grid.major.y = element_blank(),
panel.border = element_blank(),
axis.ticks.y = element_blank()) +
geom_text(aes(x = Group, y = LST24hrwMean + 0.2 * sign(LST24hrwMean),
label = paste0(round(LST24hrwMean,1),"°")),
hjust = 1.6, vjust = -.8, size = 3,
color=rgb(100,100,100, maxColorValue=255)) +
scale_y_continuous(labels = function(x) paste0(x, "°")) +
geom_hline(yintercept = LST24hrPop$LST24hrMean, linetype = "dashed") +
geom_text(aes(x = "Low Income", y = 82, label = "Above state\naverage"),
color = "gray48") +
geom_text(aes(x = "Low Income", y = 79, label = "Below state\naverage"),
color = "gray48") +
# geom_segment(aes(x = "Under 5", xend = "Under 5", y = 81, yend = 82),
# arrow = arrow(length = unit(0.3,"cm"))) +
expand_limits(y = c(78,83))
ggsave("images/MA_LST24_graph.png")
```
\pagebreak
```{r tabLST24hrPopAvgMA, fig.align = "center", fig.cap="Population-weighted average exposures to day-night average LST for priority populations in Massachusetts relative to the state average."}
LST24hrPopTab %>%
mutate(Group = recode(Group, "Minority" = "People of Color",
"MA Minority" = "MA POC",
"No HS Dip" = "No HS Diploma")) %>%
arrange(Group) %>%
kableExtra::kable(longtable = T, booktabs = T, digits = 2,
col.names = c("Group","Avg 24hr LST (°F)",
"Difference from State Avg* (°F)",
"Pct Above/Below State Avg*"),
caption = "Population-Weighted Temperature Exposure - Day-Night Average LST", align = "r") %>%
kableExtra::column_spec(2:4, width = "4cm") %>%
kableExtra::kable_styling(latex_options = c("repeat_header")) %>%
kableExtra::footnote(., symbol = paste0("State average day-night LST is ", round(mean(ma_blkgrpLST_sf$meanAvgLST,na.rm=TRUE),2), "°"))
```
```{r UHI24hrTempPop, include=FALSE}
# Pop Weighted avg of 24hr LST exposure for all Groups in Massachusetts relative to MA average
UHI24hrPop <- ma_blkgrpLST_sf %>%
as.data.frame() %>%
mutate(MA_LOWINC = if_else(MA_INCOME == "I", totalpopE, 0)) %>%
mutate(MA_LOWINC = replace_na(MA_LOWINC,0)) %>%
mutate(MA_MINORITY = if_else(MA_MINORITY == "M", totalpopE,0)) %>%
mutate(MA_MINORITY = replace_na(MA_MINORITY,0)) %>%
mutate(MA_ENGLISH = if_else(MA_ENGLISH == "E", totalpopE,0)) %>%
mutate(MA_ENGLISH = replace_na(MA_ENGLISH,0)) %>%
# filter(STATE == "Massachusetts") %>%
dplyr::select(totalpopE,
nhwhitepopE,
minorityE,
nhblackpopE,
nhamerindpopE,
nhasianpopE,
nhnativhpopE,
nhotherpopE,
nh2morepopE,
hisppopE,
povknownE,
num2povE,
eng_hhE,
eng_limitE,
age25upE,
lthsE,
allAgesE,
under5E,
over64E,
MA_LOWINC,
MA_MINORITY,
MA_ENGLISH,
UHI24avg) %>%
gather(key = Group, value = Pop, totalpopE:MA_ENGLISH) %>%
group_by(Group) %>%
summarize(UHI24hrwMean = weighted.mean(x = UHI24avg, w = Pop, na.rm = TRUE),
UHI24hrMean = mean(UHI24avg, na.rm = TRUE)) %>%
filter(Group %in% c("minorityE","eng_limitE","num2povE","lthsE","under5E","over64E","MA_LOWINC","MA_MINORITY","MA_ENGLISH")) %>%
mutate(Group = case_when(
Group == "minorityE" ~ "Minority",
Group == "eng_limitE" ~ "Limited English HH",
Group == "num2povE" ~ "Low Income",
Group == "lthsE" ~ "No HS Dip",
Group == "under5E" ~ "Under 5",
Group == "over64E" ~ "Over 64",
Group == "MA_LOWINC" ~ "MA Low Income",
Group == "MA_MINORITY" ~ "MA Minority",
Group == "MA_ENGLISH" ~ "MA Limited English HH"))
# Pop Weighted avg of PM2.5 for all Groups in Massachusetts relative to MA average
UHI24hrPopTab <- ma_blkgrpLST_sf %>%
as.data.frame() %>%
mutate(MA_LOWINC = if_else(MA_INCOME == "I", totalpopE, 0)) %>%
mutate(MA_LOWINC = replace_na(MA_LOWINC,0)) %>%
mutate(MA_MINORITY = if_else(MA_MINORITY == "M", totalpopE,0)) %>%
mutate(MA_MINORITY = replace_na(MA_MINORITY,0)) %>%
mutate(MA_ENGLISH = if_else(MA_ENGLISH == "E", totalpopE,0)) %>%
mutate(MA_ENGLISH = replace_na(MA_ENGLISH,0)) %>%
# filter(STATE == "Massachusetts") %>%
dplyr::select(totalpopE,
nhwhitepopE,
minorityE,
nhblackpopE,
nhamerindpopE,
nhasianpopE,
nhnativhpopE,
nhotherpopE,
nh2morepopE,
hisppopE,
povknownE,
num2povE,
eng_hhE,
eng_limitE,
age25upE,
lthsE,
allAgesE,
under5E,
over64E,
MA_LOWINC,
MA_MINORITY,
MA_ENGLISH,
UHI24avg) %>%
gather(key = Group, value = Pop, totalpopE:MA_ENGLISH) %>%
group_by(Group) %>%
summarize(UHI24hrwMean = weighted.mean(x = UHI24avg, w = Pop, na.rm = TRUE),
UHI24hrMean = mean(UHI24avg, na.rm = TRUE)) %>%
spread(key = Group, value = UHI24hrwMean) %>%
transmute(Minority = (minorityE/UHI24hrMean - 1)*100,
#Minority_NHW = (minorityE/nhwhitepopE - 1)*100,
`Limited English HH` = (eng_limitE/UHI24hrMean - 1)*100,
`Low Income` = (num2povE/UHI24hrMean - 1)*100,
`No HS Dip` = (lthsE/UHI24hrMean - 1)*100,
`Under 5` = (under5E/UHI24hrMean - 1)*100,
`Over 64` = (over64E/UHI24hrMean - 1)*100,