code to create annual maps of pm25 and maps of change over time

2022-02-20 23:00:32 +03:00 · 2021-12-06 12:15:46 +00:00
parent 715388dce2
commit 463bb10e3d
3 changed files with 155 additions and 0 deletions
--- a/FionaSpooner/air_pollution_pm25/annual_maps.R
+++ b/FionaSpooner/air_pollution_pm25/annual_maps.R
@@ -0,0 +1,17 @@
+library(terra)
+library(tmap)
+library(raster)
+library(RColorBrewer)
+
+source('map_func.R')
+
+#Download data from https://wustl.app.box.com/v/ACAG-V5GL02-GWRPM25/folder/148054977849
+
+years <- 1998:2020
+
+for (year in years){
+  make_pm_25(year, save_file = TRUE)
+}
+
+
+
--- a/FionaSpooner/air_pollution_pm25/change_maps.R
+++ b/FionaSpooner/air_pollution_pm25/change_maps.R
@@ -0,0 +1,103 @@
+library(terra)
+library(tmap)
+library(raster)
+
+data(World)
+world <- World
+
+pmf <- list.files('data', pattern = 'nc', full.names = TRUE)
+
+pms<- terra::rast(pmf)
+ext(pms)<- c(-180, 180, -55, 68)
+terra::crs(pms) <- "+proj=longlat +datum=WGS84 +no_defs"
+names(pms)<- 1998:2020
+
+test <- pms
+
+# Running a linear model through each cell and extracting the slope
+X <- cbind(1, 1:nlyr(test))
+## pre-computing constant part of least squares
+invXtX <- solve(t(X) %*% X) %*% t(X)
+## [2] is to just get the slope
+quickfun <- function(y) (invXtX %*% y)[2]
+#slope <- terra::app(test, quickfun) 
+
+roc <- terra::app(test, fun=quickfun, cores = 8)
+names(roc) <- "roc"
+pal = c("#1A9850","#66BD63","#A6D96A", "#D9EF8B", "#FEE08B" ,"#FDAE61" ,"#F46D43" ,"#D73027")
+
+#writeRaster(roc, "output/rate_of_change_1998_2020.tif")
+#roc <- raster::raster("output/rate_of_change_1998_2020.tif")
+
+names(roc) <- "roc"
+
+tm <- tm_shape(test2) +
+  tm_raster('roc',
+           # title = expression(atop('Average annual','change in PM'[2.5])), 
+           title = expression('Average annual change in PM'[2.5]),
+             breaks=c(-5, -2.5, -1,-0.5, 0,0.5, 1, 2.5, 5),
+            palette = pal)+
+  tm_shape(world, is.master = TRUE, bbox = c(-180,180, -56, 75))+#, c(-18000000, -7150000, 20000000,8400000)) +
+  tm_borders(lwd = 0.4,alpha = 0.5) +
+  tm_layout(main.title= paste0('Rate of change in concentration of fine particulate matter (1998-2020)'), 
+            main.title.position = c('center'), legend.title.size = 0.9, legend.text.size = 0.7)
+
+tmap_save(tm, paste0("output/rate_of_change.png"), dpi = 600)
+
+
+############### Difference between 2019 and 1998
+r98 <- raster::raster('data/V5GL02.HybridPM25.Global.199801-199812.nc') #%>% project(., "+proj=wintri")
+r19 <- raster::raster('data/V5GL02.HybridPM25.Global.201901-201912.nc')
+
+change <- r19 - r98 
+
+names(change) <- 'pm25'
+
+pal = c("#1A9850","#66BD63","#A6D96A", "#D9EF8B", "#FEE08B" ,"#FDAE61" ,"#F46D43" ,"#D73027")
+
+pm <- tm_shape(change) +
+  tm_raster('pm25',
+            title = expression('Change in PM'[2.5]), 
+            breaks=c(-125,-50, -25, 0,5, 10,20, 40, 80),
+            palette = pal)+
+  tm_shape(world, is.master = TRUE, bbox = c(-180,180, -56, 75))+#, c(-18000000, -7150000, 20000000,8400000)) +
+  tm_borders(lwd = 0.4,alpha = 0.5) +
+  tm_layout(main.title= paste0('Change in Annual Mean Concentration of Fine Particulate Matter\n 2019 vs 1998'), 
+            main.title.position = c('center'))
+
+tmap_save(pm, "Comparison_1998_2019.png")
+
+
+################ Difference between 3 yr averages around 1999 and 2019
+
+pmf <- list.files('data', pattern = 'nc', full.names = TRUE)
+early <- terra::rast(pmf[1:3])
+ext(early)<- c(-180, 180, -55, 68)
+terra::crs(early) <- "+proj=longlat +datum=WGS84 +no_defs"
+names(early)<- 1998:2000
+
+early_mean <- terra::app(early, fun=mean, cores = 8)
+
+late <- terra::rast(pmf[21:23])
+ext(late)<- c(-180, 180, -55, 68)
+terra::crs(late) <- "+proj=longlat +datum=WGS84 +no_defs"
+names(late)<- 2018:2020
+
+late_mean <- terra::app(late, fun=mean, cores = 8)
+
+difference <- late_mean - early_mean
+
+pm <- tm_shape(difference) +
+  tm_raster('mean',
+            title = expression('Change in PM'[2.5]), 
+            breaks=c(-80,-40, -20,-10,-5,0,5, 10,20, 40, 80, 100),
+            palette = pal)+
+  tm_shape(world, is.master = TRUE, bbox = c(-180,180, -56, 75))+#, c(-18000000, -7150000, 20000000,8400000)) +
+  tm_borders(lwd = 0.4,alpha = 0.5) +
+  tm_layout(main.title= paste0('Change in Annual Mean Concentration of Fine Particulate Matter\n 2019 vs 1999'), 
+            main.title.position = c('center'))
+
+pm
+
+
+
--- a/FionaSpooner/air_pollution_pm25/map_func.R
+++ b/FionaSpooner/air_pollution_pm25/map_func.R
@@ -0,0 +1,35 @@
+library(terra)
+library(tmap)
+library(raster)
+library(RColorBrewer)
+
+data(World)
+world <- World
+
+pal <- c("#8CB78A","#DDD10B","#DA970A", "#DE4F0B","#C20A08",
+          "#940B27",brewer.pal(11, "RdBu")[1],"#000000")
+
+make_pm_25 <- function(year, save_file = FALSE){
+  
+  file <- gsub("2020", year, "data/V5GL02.HybridPM25.Global.202001-202012.nc")
+  r <- raster::raster(file) 
+  names(r)<- 'pm25'
+  p <- tm_shape(r) +
+    tm_raster('pm25',
+              title = expression('PM'[2.5]), 
+              breaks=c(0, 5, 10, 15, 25, 35, 50, 100, 300),
+              palette = pal3)+
+    tm_shape(world, is.master = TRUE, bbox = c(-180,180, -56, 75))+#, c(-18000000, -7150000, 20000000,8400000)) +
+    tm_borders(lwd = 0.4,alpha = 0.5) +
+    tm_layout(main.title= paste0('Annual Mean Concentration of Fine Particulate Matter - ', year), 
+              main.title.position = c('center'))
+  if (save_file){
+    tmap_save(p, paste0("output/",year,"_global_pm25.png"))
+  }
+  
+  return(p)
+}
+
+
+
+