From 7471ea16efa1e6160c640e727b8c2785d137ca62 Mon Sep 17 00:00:00 2001 From: Tracey Mangin Date: Mon, 22 Aug 2022 12:43:35 -0700 Subject: [PATCH 01/15] test package on new zealand --- code/test-package/test_package_run.R | 30 ++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) create mode 100644 code/test-package/test_package_run.R diff --git a/code/test-package/test_package_run.R b/code/test-package/test_package_run.R new file mode 100644 index 0000000..1e217b2 --- /dev/null +++ b/code/test-package/test_package_run.R @@ -0,0 +1,30 @@ +## Tracey Mangin +## August 22, 2022 +## Test climate package + +# ## install package +# devtools::install_github("tcarleton/stagg") + +## attach libraries necessary for sourcing functions/setup +library(stagg) +library(tidyverse) +library(sf) +library(here) + + +## set paths +# root_dir <- file.path("/home/traceymangin/climate-aggregate") +input_dir <- file.path("/home/tcarleton/Climate") +save_dir <- file.path("/home/traceymangin") + +## input polygons +input_polygons <- read_sf(file.path(input_dir, "data", "shapefiles", "NZL", "gadm36_NZL_1.shp")) + +## step 1: resample a secondary weight +polygon_extent <- extent(input_polygons) +polygon_extent_vec <- c(polygon_extent@xmin, polygon_extent@xmax, polygon_extent@ymin, polygon_extent@ymax) + +cropland_weights <- secondary_weights(secondary_raster = cropland_world_2003_era5, + grid = era5_grid, + extent = polygon_extent_vec) + From 9a1f6a8ade7b10376aa3a0edb16d521f5dae0dad Mon Sep 17 00:00:00 2001 From: Tracey Mangin Date: Mon, 22 Aug 2022 12:43:48 -0700 Subject: [PATCH 02/15] test package on nz --- code/test-package/test_package_run.R | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/code/test-package/test_package_run.R b/code/test-package/test_package_run.R index 1e217b2..1aabce0 100644 --- a/code/test-package/test_package_run.R +++ b/code/test-package/test_package_run.R @@ -18,11 +18,18 @@ input_dir <- file.path("/home/tcarleton/Climate") save_dir <- file.path("/home/traceymangin") ## input polygons -input_polygons <- read_sf(file.path(input_dir, "data", "shapefiles", "NZL", "gadm36_NZL_1.shp")) + +## cluster +# input_polygons <- read_sf(file.path(input_dir, "data", "shapefiles", "NZL", "gadm36_NZL_1.shp")) + +## local +input_polygons <- read_sf('/Volumes/GoogleDrive/Shared Drives/emlab/projects/current-projects/climate-data-pipeline/shapefiles/NZL/gadm36_NZL_1.shp') + ## step 1: resample a secondary weight polygon_extent <- extent(input_polygons) polygon_extent_vec <- c(polygon_extent@xmin, polygon_extent@xmax, polygon_extent@ymin, polygon_extent@ymax) +# polygon_extent_vec <- c(-178.82695, 179.06583, -52.62083, -29.23166) cropland_weights <- secondary_weights(secondary_raster = cropland_world_2003_era5, grid = era5_grid, From a9b3a85a9f22ecf91cc34b392c8e57594baf8577 Mon Sep 17 00:00:00 2001 From: Tracey Mangin Date: Mon, 22 Aug 2022 15:09:09 -0700 Subject: [PATCH 03/15] work on testing package --- code/test-package/test_package_run.R | 43 ++++++++++++++++++++++------ 1 file changed, 35 insertions(+), 8 deletions(-) diff --git a/code/test-package/test_package_run.R b/code/test-package/test_package_run.R index 1aabce0..4f14959 100644 --- a/code/test-package/test_package_run.R +++ b/code/test-package/test_package_run.R @@ -11,27 +11,54 @@ library(tidyverse) library(sf) library(here) - ## set paths # root_dir <- file.path("/home/traceymangin/climate-aggregate") input_dir <- file.path("/home/tcarleton/Climate") save_dir <- file.path("/home/traceymangin") -## input polygons +## set inputs for run +## ---------------------------------------------------------------- + +## secondardy weights +sec_weight <- cropland_world_2003_era5 ## cropland_world_2003_era5 = crops, pop_world_2015_era5 = pop +## polygon ## cluster # input_polygons <- read_sf(file.path(input_dir, "data", "shapefiles", "NZL", "gadm36_NZL_1.shp")) ## local input_polygons <- read_sf('/Volumes/GoogleDrive/Shared Drives/emlab/projects/current-projects/climate-data-pipeline/shapefiles/NZL/gadm36_NZL_1.shp') +## polygon id +polygon_id <- 'NAME_1' + +## Step 1: filter weights for polygon extent +## ----------------------------------------------------- -## step 1: resample a secondary weight polygon_extent <- extent(input_polygons) -polygon_extent_vec <- c(polygon_extent@xmin, polygon_extent@xmax, polygon_extent@ymin, polygon_extent@ymax) -# polygon_extent_vec <- c(-178.82695, 179.06583, -52.62083, -29.23166) -cropland_weights <- secondary_weights(secondary_raster = cropland_world_2003_era5, - grid = era5_grid, - extent = polygon_extent_vec) +sec_weight_filt <- dplyr::filter(sec_weight, + x >= round(polygon_extent@xmin) - 1, + x <= round(polygon_extent@xmax) + 1, + y >= round(polygon_extent@ymin) - 1, + y <= round(polygon_extent@ymax) + 1) + +## Step 1: Overlay administrative regions onto your data's grid +## ------------------------------------------------------------ + +polygon_weights <- overlay_weights(polygons = input_polygons, + polygon_id_col = polygon_id, + grid = era5_grid, + secondary_weights = sec_weight_filt) + +## check these outputs with those from non-package run +nz_weights <- fread('/Volumes/GoogleDrive/Shared Drives/emlab/projects/current-projects/climate-data-pipeline/agg-outputs/weights/nzl_region_era5_area_crop_weights.csv') + +nz_weights <- nz_weights %>% + pivot_longer(names_to = "weight", values_to = "orig_value", w_area:weight) + +new_weights <- polygon_weights %>% + pivot_longer(names_to = "weight", values_to = "new_value", w_area:weight) %>% + left_join(nz_weights) %>% + mutate(diff = new_value - orig_value) From c104cb78aff5f1f0efb294563da456ecf69909f4 Mon Sep 17 00:00:00 2001 From: Tracey Mangin Date: Mon, 22 Aug 2022 22:53:36 -0700 Subject: [PATCH 04/15] work on testing package --- code/test-package/test_package_run.R | 71 ++++++++++++++++++++++++---- 1 file changed, 62 insertions(+), 9 deletions(-) diff --git a/code/test-package/test_package_run.R b/code/test-package/test_package_run.R index 4f14959..92a7f99 100644 --- a/code/test-package/test_package_run.R +++ b/code/test-package/test_package_run.R @@ -32,6 +32,32 @@ input_polygons <- read_sf('/Volumes/GoogleDrive/Shared Drives/emlab/projects/cur ## polygon id polygon_id <- 'NAME_1' +## years +years <- c(2009:2020) + +## prcp files +prcp_vec <- vector() + +for(i in 1:length(years)) { + + temp_name <- paste0('era5_prcp_', years[i], '.nc') + + prcp_vec[i] <- temp_name + +} + +## temp files +temp_vec <- vector() + +for(i in 1:length(years)) { + + temp_name <- paste0('era5_temp_', years[i], '.nc') + + temp_vec[i] <- temp_name + +} + + ## Step 1: filter weights for polygon extent ## ----------------------------------------------------- @@ -43,7 +69,7 @@ sec_weight_filt <- dplyr::filter(sec_weight, y >= round(polygon_extent@ymin) - 1, y <= round(polygon_extent@ymax) + 1) -## Step 1: Overlay administrative regions onto your data's grid +## Step 2: Overlay administrative regions onto your data's grid ## ------------------------------------------------------------ polygon_weights <- overlay_weights(polygons = input_polygons, @@ -51,14 +77,41 @@ polygon_weights <- overlay_weights(polygons = input_polygons, grid = era5_grid, secondary_weights = sec_weight_filt) -## check these outputs with those from non-package run -nz_weights <- fread('/Volumes/GoogleDrive/Shared Drives/emlab/projects/current-projects/climate-data-pipeline/agg-outputs/weights/nzl_region_era5_area_crop_weights.csv') +## Step 3: Aggregation (polynomial) +## ------------------------------------------------------------ + +## for cropping .nc files +min_x <- min(polygon_weights$x) - 0.5 +max_x <- max(polygon_weights$x) + 0.5 +min_y <- min(polygon_weights$y) - 0.5 +max_y <- max(polygon_weights$y) + 0.5 + +weights_ext <- raster::extent(min_x, max_x, min_y, max_y) + + + +## temperature + +temp_list <- list() + +for(i in 1:length(years)) { + + # Climate data file paths + ncpath <- file.path(input_dir, 'data/raw/temp') + nc_file <- paste0(ncpath, '/', temp_vec[i]) + + # Immediately crop to weights extent + clim_raster_tmp <- raster::crop(raster::stack(nc_file), weights_ext) + + temp_out <- staggregate_polynomial(clim_raster_tmp, + polygon_weights, + daily_agg = 'average', + time_agg = 'month', + degree = 5) + + temp_list[i] <- temp_out + +} -nz_weights <- nz_weights %>% - pivot_longer(names_to = "weight", values_to = "orig_value", w_area:weight) -new_weights <- polygon_weights %>% - pivot_longer(names_to = "weight", values_to = "new_value", w_area:weight) %>% - left_join(nz_weights) %>% - mutate(diff = new_value - orig_value) From 2ea4a8ea70713e78dde1577a818a37cd41972d17 Mon Sep 17 00:00:00 2001 From: Tracey Mangin Date: Wed, 24 Aug 2022 10:34:49 -0700 Subject: [PATCH 05/15] test package --- code/test-package/test_package_run.R | 36 +++++++++++++++++++++++++++- 1 file changed, 35 insertions(+), 1 deletion(-) diff --git a/code/test-package/test_package_run.R b/code/test-package/test_package_run.R index 92a7f99..4ae588b 100644 --- a/code/test-package/test_package_run.R +++ b/code/test-package/test_package_run.R @@ -9,13 +9,24 @@ library(stagg) library(tidyverse) library(sf) -library(here) +library(raster) +library(data.table) +library(rgdal) +# library(ncdf4) ## set paths # root_dir <- file.path("/home/traceymangin/climate-aggregate") input_dir <- file.path("/home/tcarleton/Climate") save_dir <- file.path("/home/traceymangin") +## for saving outputs +## --------------------------------------------------------------- + +data_source <- 'era5' +sec_weights <- TRUE +weight_type <- 'area_crop' +country_name <- 'new_zealand' + ## set inputs for run ## ---------------------------------------------------------------- @@ -77,6 +88,27 @@ polygon_weights <- overlay_weights(polygons = input_polygons, grid = era5_grid, secondary_weights = sec_weight_filt) +## save weight output + +weights_save_path <- file.path(save_dir, 'climate-out', 'weights') + +if(!dir.exists(weights_save_path)){ + + # If no - create it + message(crayon::yellow('Creating climate-out/weights/')) + dir.create(file.path(save_dir, "climate-out", "weights"), recursive=T) + +} + +## weights save name +weights_save_name <- paste0(paste(country_name, polygon_id, data_source, weight_type, sep="_"), ".csv") + +# Save message +message(crayon::yellow('Saving', weights_save_name, 'to', weights_save_path)) + +fwrite(polygon_weights, file = file.path(weights_save_path, weights_save_name)) + + ## Step 3: Aggregation (polynomial) ## ------------------------------------------------------------ @@ -101,6 +133,8 @@ for(i in 1:length(years)) { nc_file <- paste0(ncpath, '/', temp_vec[i]) # Immediately crop to weights extent + + # library(ncdf4) clim_raster_tmp <- raster::crop(raster::stack(nc_file), weights_ext) temp_out <- staggregate_polynomial(clim_raster_tmp, From b9a87f3ab9373eeabf86683cdf029b0f59d89ff1 Mon Sep 17 00:00:00 2001 From: traceybit Date: Wed, 24 Aug 2022 12:14:09 -0700 Subject: [PATCH 06/15] add parallelization --- code/test-package/test_package_run.R | 119 ++++++++++++++++++--------- 1 file changed, 80 insertions(+), 39 deletions(-) diff --git a/code/test-package/test_package_run.R b/code/test-package/test_package_run.R index 4ae588b..62f11e2 100644 --- a/code/test-package/test_package_run.R +++ b/code/test-package/test_package_run.R @@ -12,6 +12,7 @@ library(sf) library(raster) library(data.table) library(rgdal) +library(parallel) # library(ncdf4) ## set paths @@ -34,40 +35,13 @@ country_name <- 'new_zealand' sec_weight <- cropland_world_2003_era5 ## cropland_world_2003_era5 = crops, pop_world_2015_era5 = pop ## polygon -## cluster -# input_polygons <- read_sf(file.path(input_dir, "data", "shapefiles", "NZL", "gadm36_NZL_1.shp")) - -## local -input_polygons <- read_sf('/Volumes/GoogleDrive/Shared Drives/emlab/projects/current-projects/climate-data-pipeline/shapefiles/NZL/gadm36_NZL_1.shp') +input_polygons <- read_sf(file.path(input_dir, "data", "shapefiles", "NZL", "gadm36_NZL_1.shp")) ## polygon id polygon_id <- 'NAME_1' ## years years <- c(2009:2020) - -## prcp files -prcp_vec <- vector() - -for(i in 1:length(years)) { - - temp_name <- paste0('era5_prcp_', years[i], '.nc') - - prcp_vec[i] <- temp_name - -} - -## temp files -temp_vec <- vector() - -for(i in 1:length(years)) { - - temp_name <- paste0('era5_temp_', years[i], '.nc') - - temp_vec[i] <- temp_name - -} - ## Step 1: filter weights for polygon extent ## ----------------------------------------------------- @@ -101,7 +75,7 @@ if(!dir.exists(weights_save_path)){ } ## weights save name -weights_save_name <- paste0(paste(country_name, polygon_id, data_source, weight_type, sep="_"), ".csv") +weights_save_name <- paste0(paste(country_name, polygon_id, data_source, weight_type, sep="-"), ".csv") # Save message message(crayon::yellow('Saving', weights_save_name, 'to', weights_save_path)) @@ -120,32 +94,99 @@ max_y <- max(polygon_weights$y) + 0.5 weights_ext <- raster::extent(min_x, max_x, min_y, max_y) +## Run staggregate_polynomial function over multiple years in parallel +## -------------------------------------------------------------------------- +## Function for cropping raster and running staggregate -## temperature - -temp_list <- list() - -for(i in 1:length(years)) { - +run_stagg_year_temp <- function(year) { + # Climate data file paths ncpath <- file.path(input_dir, 'data/raw/temp') - nc_file <- paste0(ncpath, '/', temp_vec[i]) + nc_file <- paste0(ncpath, '/', 'era5_temp_', year, '.nc') # Immediately crop to weights extent - - # library(ncdf4) clim_raster_tmp <- raster::crop(raster::stack(nc_file), weights_ext) + ## run stagg for temp temp_out <- staggregate_polynomial(clim_raster_tmp, polygon_weights, daily_agg = 'average', time_agg = 'month', degree = 5) - temp_list[i] <- temp_out +} + +run_stagg_year_prcp <- function(year) { + + # Climate data file paths + ncpath <- file.path(input_dir, 'data/raw/prcp') + nc_file <- paste0(ncpath, '/', 'era5_prcp_', year, '.nc') + + # Immediately crop to weights extent + clim_raster_tmp <- raster::crop(raster::stack(nc_file), weights_ext) + + ## run stagg for prcp + prcp_out <- staggregate_polynomial(clim_raster_tmp, + polygon_weights, + daily_agg = 'sum', + time_agg = 'month', + degree = 3) } +## set up (cores, cluster) +no_cores <- parallel::detectCores() - 1 # Calculate the number of cores. Leave one in case something else needs to be done on the same computer at the same time. +cl <- makeCluster(no_cores, type = "FORK") # Initiate cluster. "FORK" means bring everything in your current environment with you. + +## run +stagg_multiyear_temp <- parLapply(cl, years, run_stagg_year_temp) +stagg_multiyear_prcp <- parLapply(cl, years, run_stagg_year_prcp) + +## stop cluster +stopCluster(cl) + +## rbind +stagg_multiyear_temp_all <- data.table::rbindlist(stagg_multiyear_temp) +stagg_multiyear_prcp_all <- data.table::rbindlist(stagg_multiyear_prcp) + +## save outputs + +output_save_path <- file.path(save_dir, "climate-out", "output") + +# Check if there is already an output folder + if(!dir.exists(output_save_path)){ + + # If no - create it + message(crayon::yellow('Creating climate-out/output')) + dir.create(file.path(save_dir, "climate-out", "output"), recursive=T) + + } + + save_name_temp <- paste0(paste(country_name, polygon_id, data_source, weight_type, + years[1], years[length(years)], 'temp', sep="-"), ".csv") + + save_name_prcp <- paste0(paste(country_name, polygon_id, data_source, weight_type, + years[1], years[length(years)], 'prcp', sep="-"), ".csv") + + + # Save message + message(crayon::yellow('Saving', save_name_temp, 'to', output_save_path)) + fwrite(stagg_multiyear_temp_all, file = file.path(output_save_path, save_name_temp)) + + message(crayon::yellow('Saving', save_name_prcp, 'to', output_save_path)) + fwrite(stagg_multiyear_prcp_all, file = file.path(output_save_path, save_name_prcp)) + + # Done + message(crayon::green('Done')) + + + + + + + + + From be2d0b1cd89371873b6af97965043639e2186ff3 Mon Sep 17 00:00:00 2001 From: traceybit Date: Wed, 24 Aug 2022 22:34:54 -0700 Subject: [PATCH 07/15] run again for chile --- code/test-package/test_package_run.R | 91 ++++++++++++++-------------- 1 file changed, 46 insertions(+), 45 deletions(-) diff --git a/code/test-package/test_package_run.R b/code/test-package/test_package_run.R index 62f11e2..d0c998d 100644 --- a/code/test-package/test_package_run.R +++ b/code/test-package/test_package_run.R @@ -13,20 +13,42 @@ library(raster) library(data.table) library(rgdal) library(parallel) -# library(ncdf4) ## set paths -# root_dir <- file.path("/home/traceymangin/climate-aggregate") input_dir <- file.path("/home/tcarleton/Climate") save_dir <- file.path("/home/traceymangin") ## for saving outputs ## --------------------------------------------------------------- +## define paths and create folders for saving weights +weights_save_path <- file.path(save_dir, 'climate-out', 'weights') + +if(!dir.exists(weights_save_path)){ + + # If no - create it + message(crayon::yellow('Creating climate-out/weights/')) + dir.create(file.path(save_dir, "climate-out", "weights"), recursive=T) + +} + +## define paths and create folders for saving outputs +output_save_path <- file.path(save_dir, "climate-out", "output") + +# Check if there is already an output folder +if(!dir.exists(output_save_path)){ + + # If no - create it + message(crayon::yellow('Creating climate-out/output')) + dir.create(file.path(save_dir, "climate-out", "output"), recursive=T) + +} + +## for naming files data_source <- 'era5' sec_weights <- TRUE weight_type <- 'area_crop' -country_name <- 'new_zealand' +country_name <- 'chile' ## set inputs for run ## ---------------------------------------------------------------- @@ -35,17 +57,16 @@ country_name <- 'new_zealand' sec_weight <- cropland_world_2003_era5 ## cropland_world_2003_era5 = crops, pop_world_2015_era5 = pop ## polygon -input_polygons <- read_sf(file.path(input_dir, "data", "shapefiles", "NZL", "gadm36_NZL_1.shp")) +input_polygons <- read_sf(file.path(input_dir, "data", "shapefiles", "CHL", "cl_comunas_geo_mortality_adm2.shp")) ## polygon id -polygon_id <- 'NAME_1' +polygon_id <- 'adm2_id' ## years -years <- c(2009:2020) +years <- c(1997:2013) ## Step 1: filter weights for polygon extent ## ----------------------------------------------------- - polygon_extent <- extent(input_polygons) sec_weight_filt <- dplyr::filter(sec_weight, @@ -56,30 +77,18 @@ sec_weight_filt <- dplyr::filter(sec_weight, ## Step 2: Overlay administrative regions onto your data's grid ## ------------------------------------------------------------ - polygon_weights <- overlay_weights(polygons = input_polygons, polygon_id_col = polygon_id, grid = era5_grid, secondary_weights = sec_weight_filt) -## save weight output - -weights_save_path <- file.path(save_dir, 'climate-out', 'weights') - -if(!dir.exists(weights_save_path)){ - - # If no - create it - message(crayon::yellow('Creating climate-out/weights/')) - dir.create(file.path(save_dir, "climate-out", "weights"), recursive=T) - -} - ## weights save name weights_save_name <- paste0(paste(country_name, polygon_id, data_source, weight_type, sep="-"), ".csv") -# Save message +## save message message(crayon::yellow('Saving', weights_save_name, 'to', weights_save_path)) +## save file fwrite(polygon_weights, file = file.path(weights_save_path, weights_save_name)) @@ -94,10 +103,10 @@ max_y <- max(polygon_weights$y) + 0.5 weights_ext <- raster::extent(min_x, max_x, min_y, max_y) -## Run staggregate_polynomial function over multiple years in parallel +## run staggregate_polynomial function over multiple years in parallel ## -------------------------------------------------------------------------- -## Function for cropping raster and running staggregate +## function for cropping raster and running staggregate run_stagg_year_temp <- function(year) { @@ -152,34 +161,26 @@ stagg_multiyear_temp_all <- data.table::rbindlist(stagg_multiyear_temp) stagg_multiyear_prcp_all <- data.table::rbindlist(stagg_multiyear_prcp) ## save outputs - -output_save_path <- file.path(save_dir, "climate-out", "output") - -# Check if there is already an output folder - if(!dir.exists(output_save_path)){ - - # If no - create it - message(crayon::yellow('Creating climate-out/output')) - dir.create(file.path(save_dir, "climate-out", "output"), recursive=T) - - } - - save_name_temp <- paste0(paste(country_name, polygon_id, data_source, weight_type, +save_name_temp <- paste0(paste(country_name, polygon_id, data_source, weight_type, years[1], years[length(years)], 'temp', sep="-"), ".csv") - save_name_prcp <- paste0(paste(country_name, polygon_id, data_source, weight_type, +save_name_prcp <- paste0(paste(country_name, polygon_id, data_source, weight_type, years[1], years[length(years)], 'prcp', sep="-"), ".csv") +## save message +message(crayon::yellow('Saving', save_name_temp, 'to', output_save_path)) + +## save output +fwrite(stagg_multiyear_temp_all, file = file.path(output_save_path, save_name_temp)) - # Save message - message(crayon::yellow('Saving', save_name_temp, 'to', output_save_path)) - fwrite(stagg_multiyear_temp_all, file = file.path(output_save_path, save_name_temp)) - - message(crayon::yellow('Saving', save_name_prcp, 'to', output_save_path)) - fwrite(stagg_multiyear_prcp_all, file = file.path(output_save_path, save_name_prcp)) +## save message +message(crayon::yellow('Saving', save_name_prcp, 'to', output_save_path)) + +## save output +fwrite(stagg_multiyear_prcp_all, file = file.path(output_save_path, save_name_prcp)) - # Done - message(crayon::green('Done')) +## fin +message(crayon::green('fin')) From b0a46f373c8093159b8857423453655ec3280320 Mon Sep 17 00:00:00 2001 From: Tracey Mangin Date: Thu, 25 Aug 2022 00:12:46 -0700 Subject: [PATCH 08/15] compare outputs - units are different --- code/test-package/compare_outputs.R | 62 ++++++++++++++++++++++++++++ code/test-package/test_package_run.R | 4 +- 2 files changed, 64 insertions(+), 2 deletions(-) create mode 100644 code/test-package/compare_outputs.R diff --git a/code/test-package/compare_outputs.R b/code/test-package/compare_outputs.R new file mode 100644 index 0000000..fe6b3cf --- /dev/null +++ b/code/test-package/compare_outputs.R @@ -0,0 +1,62 @@ +## Tracey Mangin +## August 24, 2022 +## compare outputs + +## libraries +library(tidyverse) +library(data.table) +library(rebus) + +## paths +main_path <- '/Volumes/GoogleDrive/Shared Drives/emlab/projects/current-projects/climate-data-pipeline/agg-outputs/' + +## stagg outputs +stagg_path <- paste0(main_path, 'stagg-out/') + +## read in files +weights_df <- fread(paste0(main_path, 'weights/nzl_region_era5_area_crop_weights.csv')) +weights_stagg_df <- fread(paste0(stagg_path, '/weights/new_zealand-NAME_1-era5-area_crop.csv')) + +nz_temp_df <- fread(paste0(main_path, 'output/nzl_region_era5_temp_average_2009_2020_polynomial_5_area_crop_weights.csv')) +nz_prcp_df <- fread(paste0(main_path, 'output/nzl_region_era5_prcp_sum_2009_2020_polynomial_3_area_crop_weights.csv')) + +nz_temp_stagg_df <- fread(paste0(stagg_path, 'outputs/new_zealand-NAME_1-era5-area_crop-2009-2020-temp.csv')) +nz_prcp_stagg_df <- fread(paste0(stagg_path, 'outputs/new_zealand-NAME_1-era5-area_crop-2009-2020-prcp.csv')) + + +## compare weights +## -------------------------------------------------- + +weights0 <- weights_df %>% + pivot_longer(names_to = 'type', values_to = 'value0', w_area:weight) + +weight_comp <- weights_stagg_df %>% + pivot_longer(names_to = 'type', values_to = 'value_new', w_area:weight) %>% + full_join(weights0) %>% + mutate(diff = value_new - value0) +## diff is zero + +## compare temp +## -------------------------------------------------- +temp0 <- nz_temp_df %>% + mutate(month = as.integer(str_remove(month, 'month_'))) %>% + pivot_longer(names_to = 'order', values_to = 'value0', order_1:order_5) + +temp_comp <- nz_temp_stagg_df %>% + pivot_longer(names_to = 'order', values_to = 'value_new', order_1:order_5) %>% + full_join(temp0) %>% + mutate(diff = value_new - value0) +## conversion... + +## compare prcp +## -------------------------------------------------- +prcp0 <- nz_prcp_df %>% + mutate(month = as.integer(str_remove(month, 'month_'))) %>% + pivot_longer(names_to = 'order', values_to = 'value0', order_1:order_3) + +temp_comp <- nz_temp_stagg_df %>% + pivot_longer(names_to = 'order', values_to = 'value_new', order_1:order_5) %>% + full_join(temp0) %>% + mutate(diff = value_new - value0) + + diff --git a/code/test-package/test_package_run.R b/code/test-package/test_package_run.R index d0c998d..f692bac 100644 --- a/code/test-package/test_package_run.R +++ b/code/test-package/test_package_run.R @@ -110,11 +110,11 @@ weights_ext <- raster::extent(min_x, max_x, min_y, max_y) run_stagg_year_temp <- function(year) { - # Climate data file paths + # climate data file paths ncpath <- file.path(input_dir, 'data/raw/temp') nc_file <- paste0(ncpath, '/', 'era5_temp_', year, '.nc') - # Immediately crop to weights extent + # immediately crop to weights extent clim_raster_tmp <- raster::crop(raster::stack(nc_file), weights_ext) ## run stagg for temp From d8928802105cd3a385e24804c273699d7dd18743 Mon Sep 17 00:00:00 2001 From: traceybit Date: Thu, 25 Aug 2022 12:44:35 -0700 Subject: [PATCH 09/15] convert units --- code/test-package/test_package_run.R | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/code/test-package/test_package_run.R b/code/test-package/test_package_run.R index f692bac..af972f6 100644 --- a/code/test-package/test_package_run.R +++ b/code/test-package/test_package_run.R @@ -117,6 +117,9 @@ run_stagg_year_temp <- function(year) { # immediately crop to weights extent clim_raster_tmp <- raster::crop(raster::stack(nc_file), weights_ext) + ## convert Kelvin to celcius + clim_raster_tmp <- clim_raster_tmp - 273.15 + ## run stagg for temp temp_out <- staggregate_polynomial(clim_raster_tmp, polygon_weights, @@ -128,13 +131,16 @@ run_stagg_year_temp <- function(year) { run_stagg_year_prcp <- function(year) { - # Climate data file paths + ## climate data file paths ncpath <- file.path(input_dir, 'data/raw/prcp') nc_file <- paste0(ncpath, '/', 'era5_prcp_', year, '.nc') - # Immediately crop to weights extent + ## immediately crop to weights extent clim_raster_tmp <- raster::crop(raster::stack(nc_file), weights_ext) + ## convert m to mm + clim_raster_tmp <- clim_raster_tmp * 1000 + ## run stagg for prcp prcp_out <- staggregate_polynomial(clim_raster_tmp, polygon_weights, From cfb265271a718ee83e1ef82abaf5be03003abecd Mon Sep 17 00:00:00 2001 From: traceybit Date: Thu, 25 Aug 2022 15:13:41 -0700 Subject: [PATCH 10/15] check outputs, temp seems to match but not prcp? --- code/test-package/compare_outputs.R | 22 +++++++++++----------- code/test-package/test_package_run.R | 8 ++++---- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/code/test-package/compare_outputs.R b/code/test-package/compare_outputs.R index fe6b3cf..5ccd384 100644 --- a/code/test-package/compare_outputs.R +++ b/code/test-package/compare_outputs.R @@ -14,14 +14,14 @@ main_path <- '/Volumes/GoogleDrive/Shared Drives/emlab/projects/current-projects stagg_path <- paste0(main_path, 'stagg-out/') ## read in files -weights_df <- fread(paste0(main_path, 'weights/nzl_region_era5_area_crop_weights.csv')) -weights_stagg_df <- fread(paste0(stagg_path, '/weights/new_zealand-NAME_1-era5-area_crop.csv')) +weights_df <- fread(paste0(main_path, 'weights/chl_area_era5_area_crop_weights.csv')) +weights_stagg_df <- fread(paste0(stagg_path, '/weights/chile-adm2_id-era5-area_crop.csv')) -nz_temp_df <- fread(paste0(main_path, 'output/nzl_region_era5_temp_average_2009_2020_polynomial_5_area_crop_weights.csv')) -nz_prcp_df <- fread(paste0(main_path, 'output/nzl_region_era5_prcp_sum_2009_2020_polynomial_3_area_crop_weights.csv')) +ch_temp_df <- fread(paste0(main_path, 'output/chl_area_era5_temp_average_1997_2013_polynomial_5_area_crop_weights.csv')) +ch_prcp_df <- fread(paste0(main_path, 'output/chl_area_era5_prcp_sum_1997_2013_polynomial_3_area_crop_weights.csv')) -nz_temp_stagg_df <- fread(paste0(stagg_path, 'outputs/new_zealand-NAME_1-era5-area_crop-2009-2020-temp.csv')) -nz_prcp_stagg_df <- fread(paste0(stagg_path, 'outputs/new_zealand-NAME_1-era5-area_crop-2009-2020-prcp.csv')) +ch_temp_stagg_df <- fread(paste0(stagg_path, 'outputs/chile-adm2_id-era5-area_crop-1997-2013-temp.csv')) +ch_prcp_stagg_df <- fread(paste0(stagg_path, 'outputs/chile-adm2_id-era5-area_crop-1997-2013-prcp.csv')) ## compare weights @@ -38,11 +38,11 @@ weight_comp <- weights_stagg_df %>% ## compare temp ## -------------------------------------------------- -temp0 <- nz_temp_df %>% +temp0 <- ch_temp_df %>% mutate(month = as.integer(str_remove(month, 'month_'))) %>% pivot_longer(names_to = 'order', values_to = 'value0', order_1:order_5) -temp_comp <- nz_temp_stagg_df %>% +temp_comp <- ch_temp_stagg_df %>% pivot_longer(names_to = 'order', values_to = 'value_new', order_1:order_5) %>% full_join(temp0) %>% mutate(diff = value_new - value0) @@ -50,12 +50,12 @@ temp_comp <- nz_temp_stagg_df %>% ## compare prcp ## -------------------------------------------------- -prcp0 <- nz_prcp_df %>% +prcp0 <- ch_prcp_df %>% mutate(month = as.integer(str_remove(month, 'month_'))) %>% pivot_longer(names_to = 'order', values_to = 'value0', order_1:order_3) -temp_comp <- nz_temp_stagg_df %>% - pivot_longer(names_to = 'order', values_to = 'value_new', order_1:order_5) %>% +prcp_comp <- ch_prcp_stagg_df %>% + pivot_longer(names_to = 'order', values_to = 'value_new', order_1:order_3) %>% full_join(temp0) %>% mutate(diff = value_new - value0) diff --git a/code/test-package/test_package_run.R b/code/test-package/test_package_run.R index af972f6..87c1623 100644 --- a/code/test-package/test_package_run.R +++ b/code/test-package/test_package_run.R @@ -48,7 +48,7 @@ if(!dir.exists(output_save_path)){ data_source <- 'era5' sec_weights <- TRUE weight_type <- 'area_crop' -country_name <- 'chile' +country_name <- 'new_zealand' ## set inputs for run ## ---------------------------------------------------------------- @@ -57,13 +57,13 @@ country_name <- 'chile' sec_weight <- cropland_world_2003_era5 ## cropland_world_2003_era5 = crops, pop_world_2015_era5 = pop ## polygon -input_polygons <- read_sf(file.path(input_dir, "data", "shapefiles", "CHL", "cl_comunas_geo_mortality_adm2.shp")) +input_polygons <- read_sf(file.path(input_dir, "data", "shapefiles", "NZL", "gadm36_NZL_1.shp")) ## polygon id -polygon_id <- 'adm2_id' +polygon_id <- 'NAME_1' ## years -years <- c(1997:2013) +years <- c(2009:2020) ## Step 1: filter weights for polygon extent ## ----------------------------------------------------- From 439a2d67db5d9121df3ded8ff9a50fede724c922 Mon Sep 17 00:00:00 2001 From: traceybit Date: Thu, 25 Aug 2022 22:47:26 -0700 Subject: [PATCH 11/15] add nz --- code/test-package/compare_outputs.R | 54 ++++++++++++++++++++++++++++- 1 file changed, 53 insertions(+), 1 deletion(-) diff --git a/code/test-package/compare_outputs.R b/code/test-package/compare_outputs.R index 5ccd384..a62d5fb 100644 --- a/code/test-package/compare_outputs.R +++ b/code/test-package/compare_outputs.R @@ -23,6 +23,17 @@ ch_prcp_df <- fread(paste0(main_path, 'output/chl_area_era5_prcp_sum_1997_2013_p ch_temp_stagg_df <- fread(paste0(stagg_path, 'outputs/chile-adm2_id-era5-area_crop-1997-2013-temp.csv')) ch_prcp_stagg_df <- fread(paste0(stagg_path, 'outputs/chile-adm2_id-era5-area_crop-1997-2013-prcp.csv')) +## NZL: read in files +nz_weights_df <- fread(paste0(main_path, 'weights/nzl_region_era5_area_crop_weights.csv')) +nz_weights_stagg_df <- fread(paste0(stagg_path, '/weights/new_zealand-NAME_1-era5-area_crop.csv')) + +nz_temp_df <- fread(paste0(main_path, 'output/nzl_region_era5_temp_average_2009_2020_polynomial_5_area_crop_weights.csv')) +nz_prcp_df <- fread(paste0(main_path, 'output/nzl_region_era5_prcp_sum_2009_2020_polynomial_3_area_crop_weights.csv')) + +nz_temp_stagg_df <- fread(paste0(stagg_path, 'outputs/new_zealand-NAME_1-era5-area_crop-2009-2020-temp.csv')) +nz_prcp_stagg_df <- fread(paste0(stagg_path, 'outputs/new_zealand-NAME_1-era5-area_crop-2009-2020-prcp.csv')) + + ## compare weights ## -------------------------------------------------- @@ -56,7 +67,48 @@ prcp0 <- ch_prcp_df %>% prcp_comp <- ch_prcp_stagg_df %>% pivot_longer(names_to = 'order', values_to = 'value_new', order_1:order_3) %>% - full_join(temp0) %>% + full_join(prcp0) %>% + mutate(diff = value_new - value0) + +## ---------------------------------------------------- +## NZL +## ---------------------------------------------------- + +## compare weights +## -------------------------------------------------- + +nz_weights0 <- nz_weights_df %>% + pivot_longer(names_to = 'type', values_to = 'value0', w_area:weight) + +nz_weight_comp <- nz_weights_stagg_df %>% + pivot_longer(names_to = 'type', values_to = 'value_new', w_area:weight) %>% + full_join(nz_weights0) %>% + mutate(diff = value_new - value0) +## diff is zero + +## compare temp +## -------------------------------------------------- +nz_temp0 <- nz_temp_df %>% + mutate(month = as.integer(str_remove(month, 'month_'))) %>% + pivot_longer(names_to = 'order', values_to = 'value0', order_1:order_5) + +nz_temp_comp <- nz_temp_stagg_df %>% + pivot_longer(names_to = 'order', values_to = 'value_new', order_1:order_5) %>% + full_join(nz_temp0) %>% mutate(diff = value_new - value0) +## conversion... + +## compare prcp +## -------------------------------------------------- +nz_prcp0 <- nz_prcp_df %>% + mutate(month = as.integer(str_remove(month, 'month_'))) %>% + pivot_longer(names_to = 'order', values_to = 'value0', order_1:order_3) + +nz_prcp_comp <- nz_prcp_stagg_df %>% + pivot_longer(names_to = 'order', values_to = 'value_new', order_1:order_3) %>% + full_join(nz_prcp0) %>% + mutate(diff = value_new - value0) + + From 1453a17687df3c875f49f765429aa38688239144 Mon Sep 17 00:00:00 2001 From: Tracey Mangin Date: Fri, 26 Aug 2022 00:13:20 -0700 Subject: [PATCH 12/15] add run script with stagg --- code/run_agg_stagg.R | 194 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 194 insertions(+) create mode 100644 code/run_agg_stagg.R diff --git a/code/run_agg_stagg.R b/code/run_agg_stagg.R new file mode 100644 index 0000000..6b7849f --- /dev/null +++ b/code/run_agg_stagg.R @@ -0,0 +1,194 @@ +## Tracey Mangin +## August 25, 2022 +## Workflow using stagg package + +# ## install package +# devtools::install_github("tcarleton/stagg") + +## attach libraries necessary for sourcing functions/setup +library(stagg) +library(tidyverse) +library(sf) +library(raster) +library(data.table) +library(rgdal) +library(parallel) + +## set paths +input_dir <- file.path("/home/tcarleton/Climate") ## path for shapefiles and raw climate data +save_dir <- file.path("/home/traceymangin") ## path for saving outputs (folders defined below) + +## for saving outputs +## --------------------------------------------------------------- + +## define paths and create folders for saving weights +weights_save_path <- file.path(save_dir, 'climate-out', 'weights') + +if(!dir.exists(weights_save_path)){ + + # If no - create it + message(crayon::yellow('Creating climate-out/weights/')) + dir.create(file.path(save_dir, "climate-out", "weights"), recursive=T) + +} + +## define paths and create folders for saving outputs +output_save_path <- file.path(save_dir, "climate-out", "output") + +# Check if there is already an output folder +if(!dir.exists(output_save_path)){ + + # If no - create it + message(crayon::yellow('Creating climate-out/output')) + dir.create(file.path(save_dir, "climate-out", "output"), recursive=T) + +} + +## for naming files +data_source <- 'era5' +sec_weights <- TRUE +weight_type <- 'area_crop' +country_name <- 'new_zealand' + +## set inputs for run +## ---------------------------------------------------------------- + +## secondary weights (data included in the stagg package, no need to load separately) +## cropland_world_2003_era5 = crops +## pop_world_2015_era5 = pop +sec_weight <- cropland_world_2003_era5 + +## polygon +input_polygons <- read_sf(file.path(input_dir, "data", "shapefiles", "NZL", "gadm36_NZL_1.shp")) + +## polygon id +polygon_id <- 'NAME_1' + +## years +years <- c(2009:2020) + +## Step 1: filter weights for polygon extent +## ----------------------------------------------------- +polygon_extent <- extent(input_polygons) + +sec_weight_filt <- dplyr::filter(sec_weight, + x >= round(polygon_extent@xmin) - 1, + x <= round(polygon_extent@xmax) + 1, + y >= round(polygon_extent@ymin) - 1, + y <= round(polygon_extent@ymax) + 1) + +## Step 2: Overlay administrative regions onto your data's grid +## ------------------------------------------------------------ +polygon_weights <- overlay_weights(polygons = input_polygons, + polygon_id_col = polygon_id, + grid = era5_grid, + secondary_weights = sec_weight_filt) + +## weights save name +weights_save_name <- paste0(paste(country_name, polygon_id, data_source, weight_type, sep="-"), ".csv") + +## save message +message(crayon::yellow('Saving', weights_save_name, 'to', weights_save_path)) + +## save file +fwrite(polygon_weights, file = file.path(weights_save_path, weights_save_name)) + + +## Step 3: Aggregation (polynomial) +## ------------------------------------------------------------ + +## for cropping .nc files +min_x <- min(polygon_weights$x) - 0.5 +max_x <- max(polygon_weights$x) + 0.5 +min_y <- min(polygon_weights$y) - 0.5 +max_y <- max(polygon_weights$y) + 0.5 + +weights_ext <- raster::extent(min_x, max_x, min_y, max_y) + +## run staggregate_polynomial function over multiple years in parallel +## -------------------------------------------------------------------------- + +## function for cropping raster and running staggregate + +run_stagg_year_temp <- function(year) { + + # climate data file paths + ncpath <- file.path(input_dir, 'data/raw/temp') + nc_file <- paste0(ncpath, '/', 'era5_temp_', year, '.nc') + + # immediately crop to weights extent + clim_raster_tmp <- raster::crop(raster::stack(nc_file), weights_ext) + + ## convert Kelvin to celcius + clim_raster_tmp <- clim_raster_tmp - 273.15 + + ## run stagg for temp + temp_out <- staggregate_polynomial(clim_raster_tmp, + polygon_weights, + daily_agg = 'average', + time_agg = 'month', + degree = 5) + +} + +run_stagg_year_prcp <- function(year) { + + ## climate data file paths + ncpath <- file.path(input_dir, 'data/raw/prcp') + nc_file <- paste0(ncpath, '/', 'era5_prcp_', year, '.nc') + + ## immediately crop to weights extent + clim_raster_tmp <- raster::crop(raster::stack(nc_file), weights_ext) + + ## convert m to mm + clim_raster_tmp <- clim_raster_tmp * 1000 + + ## run stagg for prcp + prcp_out <- staggregate_polynomial(clim_raster_tmp, + polygon_weights, + daily_agg = 'sum', + time_agg = 'month', + degree = 3) + +} + + +## set up (cores, cluster) +no_cores <- parallel::detectCores() - 1 # Calculate the number of cores. Leave one in case something else needs to be done on the same computer at the same time. +cl <- makeCluster(no_cores, type = "FORK") # Initiate cluster. "FORK" means bring everything in your current environment with you. + +## run +stagg_multiyear_temp <- parLapply(cl, years, run_stagg_year_temp) +stagg_multiyear_prcp <- parLapply(cl, years, run_stagg_year_prcp) + +## stop cluster +stopCluster(cl) + +## rbind +stagg_multiyear_temp_all <- data.table::rbindlist(stagg_multiyear_temp) +stagg_multiyear_prcp_all <- data.table::rbindlist(stagg_multiyear_prcp) + +## save outputs +save_name_temp <- paste0(paste(country_name, polygon_id, data_source, weight_type, + years[1], years[length(years)], 'temp', sep="-"), ".csv") + +save_name_prcp <- paste0(paste(country_name, polygon_id, data_source, weight_type, + years[1], years[length(years)], 'prcp', sep="-"), ".csv") + +## save message +message(crayon::yellow('Saving', save_name_temp, 'to', output_save_path)) + +## save output +fwrite(stagg_multiyear_temp_all, file = file.path(output_save_path, save_name_temp)) + +## save message +message(crayon::yellow('Saving', save_name_prcp, 'to', output_save_path)) + +## save output +fwrite(stagg_multiyear_prcp_all, file = file.path(output_save_path, save_name_prcp)) + +## fin +message(crayon::green('fin')) + + + From d0ea521e4e92a750d4d00e5120662a17bdb54adb Mon Sep 17 00:00:00 2001 From: traceybit Date: Fri, 26 Aug 2022 10:46:17 -0700 Subject: [PATCH 13/15] adjust to use input df --- code/run_agg_stagg.R | 40 ++++++++++++++++++++++++++-------------- 1 file changed, 26 insertions(+), 14 deletions(-) diff --git a/code/run_agg_stagg.R b/code/run_agg_stagg.R index 6b7849f..5212073 100644 --- a/code/run_agg_stagg.R +++ b/code/run_agg_stagg.R @@ -16,6 +16,7 @@ library(parallel) ## set paths input_dir <- file.path("/home/tcarleton/Climate") ## path for shapefiles and raw climate data +country_input_dir <- file.path("/home/traceymangin/data/inputs") ## path for country inputs save_dir <- file.path("/home/traceymangin") ## path for saving outputs (folders defined below) ## for saving outputs @@ -44,28 +45,39 @@ if(!dir.exists(output_save_path)){ } -## for naming files -data_source <- 'era5' -sec_weights <- TRUE -weight_type <- 'area_crop' -country_name <- 'new_zealand' -## set inputs for run -## ---------------------------------------------------------------- +## read in inputs and define variables for filtering/naming/running +## ----------------------------------------------------------- ## secondary weights (data included in the stagg package, no need to load separately) ## cropland_world_2003_era5 = crops ## pop_world_2015_era5 = pop sec_weight <- cropland_world_2003_era5 -## polygon -input_polygons <- read_sf(file.path(input_dir, "data", "shapefiles", "NZL", "gadm36_NZL_1.shp")) +## for naming files +weight_type <- 'area_crop' + +## read country inputs and filter for country that you want to run +country_inputs <- fread(file.path(country_input_dir, "climate_country_inputs.csv")) -## polygon id -polygon_id <- 'NAME_1' +## choose country that you want to run +country_name <- 'ECU' ## use country abbreviation (see country column in country_inputs) + +## filter inputs for country, define items +data_source <- country_inputs[country == country_name, data_src][1] + +## read in polygon +poly_name <- country_inputs[country == country_name, shapefile_name][1] +input_polygons <- read_sf(file.path(input_dir, "data", "shapefiles", country_name, poly_name)) + +## polygon geoid +polygon_id <- country_inputs[country == country_name, id_var][1] ## years -years <- c(2009:2020) +min_year <- country_inputs[country == country_name, start_year][1] +max_year <- country_inputs[country == country_name, end_year][1] +years <- c(min_year:max_year) + ## Step 1: filter weights for polygon extent ## ----------------------------------------------------- @@ -170,10 +182,10 @@ stagg_multiyear_prcp_all <- data.table::rbindlist(stagg_multiyear_prcp) ## save outputs save_name_temp <- paste0(paste(country_name, polygon_id, data_source, weight_type, - years[1], years[length(years)], 'temp', sep="-"), ".csv") + min_year, max_year, 'temp', sep="-"), ".csv") save_name_prcp <- paste0(paste(country_name, polygon_id, data_source, weight_type, - years[1], years[length(years)], 'prcp', sep="-"), ".csv") + min_year, max_year, 'prcp', sep="-"), ".csv") ## save message message(crayon::yellow('Saving', save_name_temp, 'to', output_save_path)) From 9fcbcb309ef6a45221da8eae0cbe1b45c145a7bf Mon Sep 17 00:00:00 2001 From: traceybit Date: Fri, 26 Aug 2022 13:13:13 -0700 Subject: [PATCH 14/15] delete trial script --- code/test-package/test_package_run.R | 199 --------------------------- 1 file changed, 199 deletions(-) delete mode 100644 code/test-package/test_package_run.R diff --git a/code/test-package/test_package_run.R b/code/test-package/test_package_run.R deleted file mode 100644 index 87c1623..0000000 --- a/code/test-package/test_package_run.R +++ /dev/null @@ -1,199 +0,0 @@ -## Tracey Mangin -## August 22, 2022 -## Test climate package - -# ## install package -# devtools::install_github("tcarleton/stagg") - -## attach libraries necessary for sourcing functions/setup -library(stagg) -library(tidyverse) -library(sf) -library(raster) -library(data.table) -library(rgdal) -library(parallel) - -## set paths -input_dir <- file.path("/home/tcarleton/Climate") -save_dir <- file.path("/home/traceymangin") - -## for saving outputs -## --------------------------------------------------------------- - -## define paths and create folders for saving weights -weights_save_path <- file.path(save_dir, 'climate-out', 'weights') - -if(!dir.exists(weights_save_path)){ - - # If no - create it - message(crayon::yellow('Creating climate-out/weights/')) - dir.create(file.path(save_dir, "climate-out", "weights"), recursive=T) - -} - -## define paths and create folders for saving outputs -output_save_path <- file.path(save_dir, "climate-out", "output") - -# Check if there is already an output folder -if(!dir.exists(output_save_path)){ - - # If no - create it - message(crayon::yellow('Creating climate-out/output')) - dir.create(file.path(save_dir, "climate-out", "output"), recursive=T) - -} - -## for naming files -data_source <- 'era5' -sec_weights <- TRUE -weight_type <- 'area_crop' -country_name <- 'new_zealand' - -## set inputs for run -## ---------------------------------------------------------------- - -## secondardy weights -sec_weight <- cropland_world_2003_era5 ## cropland_world_2003_era5 = crops, pop_world_2015_era5 = pop - -## polygon -input_polygons <- read_sf(file.path(input_dir, "data", "shapefiles", "NZL", "gadm36_NZL_1.shp")) - -## polygon id -polygon_id <- 'NAME_1' - -## years -years <- c(2009:2020) - -## Step 1: filter weights for polygon extent -## ----------------------------------------------------- -polygon_extent <- extent(input_polygons) - -sec_weight_filt <- dplyr::filter(sec_weight, - x >= round(polygon_extent@xmin) - 1, - x <= round(polygon_extent@xmax) + 1, - y >= round(polygon_extent@ymin) - 1, - y <= round(polygon_extent@ymax) + 1) - -## Step 2: Overlay administrative regions onto your data's grid -## ------------------------------------------------------------ -polygon_weights <- overlay_weights(polygons = input_polygons, - polygon_id_col = polygon_id, - grid = era5_grid, - secondary_weights = sec_weight_filt) - -## weights save name -weights_save_name <- paste0(paste(country_name, polygon_id, data_source, weight_type, sep="-"), ".csv") - -## save message -message(crayon::yellow('Saving', weights_save_name, 'to', weights_save_path)) - -## save file -fwrite(polygon_weights, file = file.path(weights_save_path, weights_save_name)) - - -## Step 3: Aggregation (polynomial) -## ------------------------------------------------------------ - -## for cropping .nc files -min_x <- min(polygon_weights$x) - 0.5 -max_x <- max(polygon_weights$x) + 0.5 -min_y <- min(polygon_weights$y) - 0.5 -max_y <- max(polygon_weights$y) + 0.5 - -weights_ext <- raster::extent(min_x, max_x, min_y, max_y) - -## run staggregate_polynomial function over multiple years in parallel -## -------------------------------------------------------------------------- - -## function for cropping raster and running staggregate - -run_stagg_year_temp <- function(year) { - - # climate data file paths - ncpath <- file.path(input_dir, 'data/raw/temp') - nc_file <- paste0(ncpath, '/', 'era5_temp_', year, '.nc') - - # immediately crop to weights extent - clim_raster_tmp <- raster::crop(raster::stack(nc_file), weights_ext) - - ## convert Kelvin to celcius - clim_raster_tmp <- clim_raster_tmp - 273.15 - - ## run stagg for temp - temp_out <- staggregate_polynomial(clim_raster_tmp, - polygon_weights, - daily_agg = 'average', - time_agg = 'month', - degree = 5) - -} - -run_stagg_year_prcp <- function(year) { - - ## climate data file paths - ncpath <- file.path(input_dir, 'data/raw/prcp') - nc_file <- paste0(ncpath, '/', 'era5_prcp_', year, '.nc') - - ## immediately crop to weights extent - clim_raster_tmp <- raster::crop(raster::stack(nc_file), weights_ext) - - ## convert m to mm - clim_raster_tmp <- clim_raster_tmp * 1000 - - ## run stagg for prcp - prcp_out <- staggregate_polynomial(clim_raster_tmp, - polygon_weights, - daily_agg = 'sum', - time_agg = 'month', - degree = 3) - -} - - -## set up (cores, cluster) -no_cores <- parallel::detectCores() - 1 # Calculate the number of cores. Leave one in case something else needs to be done on the same computer at the same time. -cl <- makeCluster(no_cores, type = "FORK") # Initiate cluster. "FORK" means bring everything in your current environment with you. - -## run -stagg_multiyear_temp <- parLapply(cl, years, run_stagg_year_temp) -stagg_multiyear_prcp <- parLapply(cl, years, run_stagg_year_prcp) - -## stop cluster -stopCluster(cl) - -## rbind -stagg_multiyear_temp_all <- data.table::rbindlist(stagg_multiyear_temp) -stagg_multiyear_prcp_all <- data.table::rbindlist(stagg_multiyear_prcp) - -## save outputs -save_name_temp <- paste0(paste(country_name, polygon_id, data_source, weight_type, - years[1], years[length(years)], 'temp', sep="-"), ".csv") - -save_name_prcp <- paste0(paste(country_name, polygon_id, data_source, weight_type, - years[1], years[length(years)], 'prcp', sep="-"), ".csv") - -## save message -message(crayon::yellow('Saving', save_name_temp, 'to', output_save_path)) - -## save output -fwrite(stagg_multiyear_temp_all, file = file.path(output_save_path, save_name_temp)) - -## save message -message(crayon::yellow('Saving', save_name_prcp, 'to', output_save_path)) - -## save output -fwrite(stagg_multiyear_prcp_all, file = file.path(output_save_path, save_name_prcp)) - -## fin -message(crayon::green('fin')) - - - - - - - - - - From 78cd08153c611ddfb60bb25dc3845bf12975deb3 Mon Sep 17 00:00:00 2001 From: traceybit Date: Fri, 26 Aug 2022 13:13:43 -0700 Subject: [PATCH 15/15] compare outputs --- code/test-package/compare_outputs.R | 46 +++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/code/test-package/compare_outputs.R b/code/test-package/compare_outputs.R index a62d5fb..bb9de91 100644 --- a/code/test-package/compare_outputs.R +++ b/code/test-package/compare_outputs.R @@ -34,6 +34,15 @@ nz_temp_stagg_df <- fread(paste0(stagg_path, 'outputs/new_zealand-NAME_1-era5-ar nz_prcp_stagg_df <- fread(paste0(stagg_path, 'outputs/new_zealand-NAME_1-era5-area_crop-2009-2020-prcp.csv')) +## ECU: read in files +ecu_weights_df <- fread(paste0(main_path, 'weights/ecu_province_era5_area_crop_weights.csv')) +ecu_weights_stagg_df <- fread(paste0(stagg_path, '/weights/ECU-NAME_1-era5-area_crop.csv')) + +ecu_temp_df <- fread(paste0(main_path, 'output/ecu_province_era5_temp_average_1990_2018_polynomial_5_area_crop_weights.csv')) +ecu_prcp_df <- fread(paste0(main_path, 'output/ecu_province_era5_prcp_sum_1990_2018_polynomial_3_area_crop_weights.csv')) + +ecu_temp_stagg_df <- fread(paste0(stagg_path, 'outputs/ECU-NAME_1-era5-area_crop-1990-2018-temp.csv')) +ecu_prcp_stagg_df <- fread(paste0(stagg_path, 'outputs/ECU-NAME_1-era5-area_crop-1990-2018-prcp.csv')) ## compare weights ## -------------------------------------------------- @@ -110,5 +119,42 @@ nz_prcp_comp <- nz_prcp_stagg_df %>% mutate(diff = value_new - value0) +## ---------------------------------------------------- +## ECU +## ---------------------------------------------------- + +## compare weights +## -------------------------------------------------- +ecu_weights0 <- ecu_weights_df %>% + pivot_longer(names_to = 'type', values_to = 'value0', w_area:weight) + +ecu_weight_comp <- ecu_weights_stagg_df %>% + pivot_longer(names_to = 'type', values_to = 'value_new', w_area:weight) %>% + full_join(ecu_weights0) %>% + mutate(diff = value_new - value0) +## diff is zero + +## compare temp +## -------------------------------------------------- +ecu_temp0 <- ecu_temp_df %>% + mutate(month = as.integer(str_remove(month, 'month_'))) %>% + pivot_longer(names_to = 'order', values_to = 'value0', order_1:order_5) + +ecu_temp_comp <- ecu_temp_stagg_df %>% + pivot_longer(names_to = 'order', values_to = 'value_new', order_1:order_5) %>% + full_join(ecu_temp0) %>% + mutate(diff = value_new - value0) +## conversion... + +## compare prcp +## -------------------------------------------------- +ecu_prcp0 <- ecu_prcp_df %>% + mutate(month = as.integer(str_remove(month, 'month_'))) %>% + pivot_longer(names_to = 'order', values_to = 'value0', order_1:order_3) + +ecu_prcp_comp <- ecu_prcp_stagg_df %>% + pivot_longer(names_to = 'order', values_to = 'value_new', order_1:order_3) %>% + full_join(ecu_prcp0) %>% + mutate(diff = value_new - value0)