diff --git a/R/map_covariates.R b/R/map_covariates.R deleted file mode 100755 index e69de29..0000000 diff --git a/R/process_spatial_covariates.R b/R/process_spatial_covariates.R index 2bc0386..d6a9ea9 100755 --- a/R/process_spatial_covariates.R +++ b/R/process_spatial_covariates.R @@ -1,35 +1,25 @@ -#' Project and save imperviousness data to a new crs -#' (used for projection to EPSG:5070 previously to prediction grid creation +#' Aggregate raster and store it in a new file (useful for dem) #' -#' @param crs_dest A character with the crs of destination -#' @param crs_name A character for the crs without weird characters -#' to concat to file name -project_imp <- function(crs_dest, crs_name) { - # imperviousness is by default in wgs84 (epsg:4326) - imp <- terra::rast("../input/NC_imperviousness_2019.tif") - if (!terra::same.crs(terra::crs(imp), crs_dest)) { - imp_proj <- terra::project(imp, crs_dest) - terra::writeRaster(imp_proj, - filename = paste0( - "../input/", - "NC_imperviousness_2019_", - crs_name, - ".tif" - ), +#' @param in_filepath A character path to data file with high resolution +#' @param out_filepath A character path to the folder where the aggregate file +#' will be stored +agg_rast <- function(in_filepath, out_filepath, agg_fact = 30) { + raw <- terra::rast(in_filepath) + agg <- terra::aggregate(raw, fact = agg_fact, fun = "median") %>% + terra::writeRaster( + filename = paste0(out_filepath), overwrite = TRUE ) - } + return(agg) } -#' Aggregate digital elevation model (~1m to ~30m) + +#' Subset a polygon area from a SpatRaster or a SpatVector #' -#' @param dem_path A path to dem file with high resolution -agg_dem <- function(dem_path = "../input/NC-DEM.tif") { - dem <- terra::rast(dem_path) - dem_agg <- terra::aggregate(dem, fact = 30, fun = "median") - terra::writeRaster(dem_agg, - filename = "../input/NC-DEM-agg.tif", - overwrite = TRUE - ) - return(dem_agg) +#' @param sp a SpatRaster or a SpatVector +#' @param poly a SpatVector with polygon geometry +subset_area <- function(sp, poly) { + poly_proj <- terra::project(poly, terra::crs(sp)) + crop_sp <- terra::crop(sp, poly_proj) + return(crop_sp) } diff --git a/R/create_prediction_grid.R b/archive_source/create_prediction_grid.R similarity index 100% rename from R/create_prediction_grid.R rename to archive_source/create_prediction_grid.R diff --git a/tests/testthat/test-create_prediction_grid.R b/tests/testthat/test-create_prediction_grid.R deleted file mode 100755 index e69de29..0000000 diff --git a/tests/testthat/test-manipulate_spacetime_data.R b/tests/testthat/test-manipulate_spacetime_data.R deleted file mode 100755 index e69de29..0000000 diff --git a/tests/testthat/test-process_monitors.R b/tests/testthat/test-process_monitors.R deleted file mode 100755 index e69de29..0000000 diff --git a/vignettes/add_covariates_to_prediction_grid.Rmd b/vignettes/add_covariates_to_prediction_grid.Rmd deleted file mode 100755 index 84da6e6..0000000 --- a/vignettes/add_covariates_to_prediction_grid.Rmd +++ /dev/null @@ -1,272 +0,0 @@ ---- -title: "Add covariates: example on RTP area" -output: html_document -date: "2023-11-09" ---- - -Load tools - -```{r, message=F} -pkgs <- c( - "tidyverse", - "terra", - "sf", - "sftime", - "reshape2", - "stringr", - "dplyr", - "exactextractr", - "tidyterra", - "viridis", - "ggplot2" -) -sapply(pkgs, library, character.only = TRUE) -source("../R/add_covariates.R") -source("../R/manipulate_spacetime_data.R") -``` - - -```{r} -nc <- vect(paste0( - "../input/NC_county_boundary/", - "North_Carolina_State_and_County_Boundary_Polygons.shp" -)) -same.crs(crs(nc), "EPSG:6543") -``` - -# Create prediction grid based on imperviousness aggregation - -```{r} -imp <- rast("../input/NC_imperviousness_2019.tif") -nc_proj <- terra::project(nc, crs(imp)) -imp <- terra::mask(imp, nc_proj) -# aggregation at 300m -pred_rast <- aggregate(imp, fact = 10, fun = "mean") -names(pred_rast) <- "imp" -plot(pred_rast) -``` -# Covariate extraction on RTP area - -### Zoom on RTP area - -```{r} -lat <- c(35.6, 36.11, 36.11, 35.6) -lon <- c(-79.19, -79.19, -78.39, -78.39) -ext_rtp <- vect(cbind(lon, lat), type = "polygons", crs = "EPSG:4326") -ext_rtp <- terra::project(ext_rtp, crs(pred_rast)) -pred_rast_rtp <- terra::crop(pred_rast, ext_rtp) -pred_vect_rtp <- terra::as.points(pred_rast_rtp) -``` - - -### Add spatial covariates - -Tree canopy height - -```{r} -pred_vect_rtp <- add_canopy_h(sp_vect = pred_vect_rtp) -``` - -Digital elevation model - -```{r} -pred_vect_rtp <- add_dem(sp_vect = pred_vect_rtp) -``` - -Tree canopy cover - -```{r} -pred_vect_rtp <- add_tcc(sp_vect = pred_vect_rtp) -``` - -Building footprint - -```{r} -pred_vect_rtp <- add_build_fp(sp_vect = pred_vect_rtp) -``` - -Building height - -```{r} -pred_vect_rtp <- add_build_h(sp_vect = pred_vect_rtp) -``` - -Land cover ratio - -```{r} -lc <- rast("../input/NC_nlcd_crs-wgs84.tif") -pred_vect_rtp <- add_nlcd_class_ratio(pred_vect_rtp, nlcd = lc) -``` - - -```{r} -pred_rast_rtp <- terra::rasterize(pred_vect_rtp, - pred_rast_rtp, - field = names(pred_vect_rtp) -) -head(pred_rast_rtp) -plot(pred_rast_rtp[[1:6]]) -plot(pred_rast_rtp[[7:21]]) -``` - -Add terrain covariates to raster - -```{r} -pred_rast_rtp <- add_terrain(pred_rast_rtp) -plot(pred_rast_rtp[[22:25]]) -``` - -### Create a SpatRasterDataset with all spatial covariates - -```{r} -pred_rastds <- list() -for (i in names(pred_rast_rtp)) { - pred_rastds[[i]] <- pred_rast_rtp[[i]] -} -pred_rastds <- terra::sds(pred_rastds) -``` - -### Add temporal covariates - -```{r} -era5 <- fread("../input/era5_daily_reanalysis_2022-05-02_2022-09-29.csv") -era5 <- era5 %>% rename(time = date) - -# convert to stdt and then to SpatRasterDataset -era5_stdt <- create_stdtobj(era5, "EPSG:4326") -era5_rastds <- convert_stdt_spatrastdataset(era5_stdt) -plot(era5_rastds$TNwmo$`2022-05-02`) - -# empty prediction SpatVector -new_pred_vect_rtp <- vect(geom(pred_vect_rtp)[, c("x", "y")], - type = "points", - crs = crs(pred_vect_rtp) -) -pred_rds_era5 <- list() -for (i in 2:7) { - pred_rds_era5[[i - 1]] <- terra::project( - new_pred_vect_rtp, - crs(era5_rastds[[i]]) - ) %>% - terra::extract(x = era5_rastds[[i]], y = ., bind = TRUE) %>% - terra::project(., crs(new_pred_vect_rtp)) %>% - terra::rasterize(., pred_rast_rtp, field = names(.)) -} -pred_rds_era5 <- terra::sds(pred_rds_era5) -names(pred_rds_era5) <- names(era5_rastds[[2:7]]) -plot(pred_rds_era5$TNwmo$`2022-05-07`) -``` -### TODEL - -```{r} -# empty prediction SpatVector -new_pred_vect_rtp <- vect(geom(pred_vect_rtp)[, c("x", "y")], - type = "points", - crs = crs(pred_vect_rtp) -) -pred_rds_era5 <- list() -for (i in 2:7) { - pred_rds_era5[[i - 1]] <- terra::project( - new_pred_vect_rtp, - crs(era5_rastds[[i]]) - ) %>% - terra::extract(x = era5_rastds[[i]], y = ., bind = TRUE) %>% - terra::project(., crs(new_pred_vect_rtp)) %>% - as.data.frame(., geom = "XY") %>% - reshape2::melt(., - id.vars = c("x", "y"), - variable.name = "time", - value.name = names(era5_rastds)[[i]] - ) -} - -df <- pred_vect_rtp %>% - as.data.frame(., geom = "XY") - - -# PROBLEME : les coordonnees ne semblent plus correspondre, peut etre a cause -# des reprojections successives -output_vect <- pred_rds_era5 %>% - reduce(inner_join, by = c("x", "y", "time")) %>% - reduce(.x = list(., df), .f = inner_join, by = c("x", "y")) %>% - rename("lon" = "x") %>% - rename("lat" = "y") %>% - vect(., - geom = c("lon", "lat"), - crs = crs(new_pred_vect_rtp), - keepgeom = TRUE - ) - - -v <- df[which(unique(hello[, c("x", "y")]) == unique(df[, c("x", "y")])), ] -nrow(v) -plot(v$x, v$y, cex = 0.01) -``` - - -### Merge both SpatRasterDatasets to have all the covariates in the same one - -```{r} -pred_rastds <- c(pred_rds_era5, pred_rastds) -plot(pred_rastds$imp) -plot(pred_rastds$frac_42) -plot(pred_rastds$TNwmo$`2022-05-27`) -plot(pred_rastds$TNwmo$`2022-08-27`) -``` - -# Covarites on entire North Carolina - -Output of create_prediction_grid_full.R (launched on geo HPC on Nov 22nd 2023) - -```{r} -eg1 <- terra::rast("../input/pred_grid_2022-05-18.tif") -eg2 <- terra::rast("../input/pred_grid_2022-08-18.tif") -``` - -Land cover index - -```{r} -nlcd_classes <- list( - value = c( - 11, 21, 22, 23, 24, 31, 41, 42, 43, 52, - 71, 81, 82, 90, 95 - ), - class = c( - "WTR", "OSD", "LID", "MID", "HID", - "BRN", "DFO", "EFO", "MFO", "SHB", - "GRS", "PAS", "CRP", "WDW", "EHW" - ), - names = c( - "Open Water", - "Developed, Open Space", - "Developed, Low Intensity", - "Developed, Medium Intensity", - "Developed, High Intensity", - "Barren Land", - "Deciduous Forest", - "Evergreen Forest", - "Mixed Forest", - "Shrub/Scrub", - "Herbaceous", - "Hay/Pasture", - "Cultivated Crops", - "Woody Wetlands", - "Emergent Herbaceous Wetlands" - ), - col = c( - "#476ba1", "#decaca", "#d99482", "#ee0000", - "#ab0000", "#b3aea3", "#68ab63", "#1c6330", - "#b5ca8f", "#ccba7d", "#e3e3c2", "#dcd93d", - "#ab7028", "#bad9eb", "#70a3ba" - ) -) -nlcd_classes <- as.data.frame(nlcd_classes) -``` - -```{r} -terra::plot(eg1[[c(1:6, 22:25)]]) -terra::plot(eg1[[7:21]]) -nlcd_classes -terra::plot(eg1[[26:31]]) -terra::plot(eg2[[26:31]]) -``` diff --git a/vignettes/create-prediction-grid.Rmd b/vignettes/create-prediction-grid.Rmd deleted file mode 100755 index 83c4bbc..0000000 --- a/vignettes/create-prediction-grid.Rmd +++ /dev/null @@ -1,470 +0,0 @@ ---- -title: "create-prediction-grid" -output: - pdf_document: default - html_document: default -date: "2023-07-20" ---- - -Libraries - -```{r, message=F} -library(terra) -library(sf) -library(ggplot2) -library(ggspatial) -library(tidyterra) -library(rgeos) -library(data.table) # -- for large flat datasets -``` - -### In this script... - -#### 1) Create empty grid - -1. Open NC borders and 30m-imperviousness data over NC -2. Compute imperviousness mean over 5km-resolution grid -3. Create an urban mask defined as 5kmx5km grid cells with imperviousness \>5% -4. Create a grid with 1kmx1km cells above rural areas and 200mx200m cells above urban ones - -#### 2) Add covariates - -- Land cover and urban indicators: imperviousness, building footprint, building height, tree canopy cover -- Digital Elevation Model (DEM) -- Meteorological covariates - -#### 3) Save prediction grid with space-time covariates in a datatable - -```{r} -input_path <- "../input/" -``` - -## Empty grid creation - -### NC borders - -```{r} -nc_poly_file <- paste0("NC_county_boundary/", - "North_Carolina_State_and_County_Boundary_Polygons.shp") -nc_borders <- vect(input_path, nc_poly_file) -``` - -Change nc_borders units (us feet to meter) - -```{r} -nc_borders <- project(nc_borders, "EPSG:5070") -``` - -### 5km-resolution imperviousness grid - -Open and map imperviousness 30m-resolution data across NC - -```{r} -imp <- rast("../input/NC_imperviousness_2019.tif") -imp <- project(imp, "EPSG:5070") - -imp_plot <- ifel(imp == 0, NA, imp) -ggplot() + - geom_spatraster(data = imp_plot) + - geom_sf( - data = st_as_sf(nc_borders), aes(geometry = geometry), - colour = "grey", linewidth = .3, fill = NA - ) + - scale_fill_whitebox_c( - palette = "muted", - labels = scales::label_number(suffix = "%"), - n.breaks = 12, - guide = guide_legend(reverse = TRUE), - na.value = NA - ) + - labs( - fill = "", - title = "Prediction grid" - ) + - annotation_scale( - location = "bl", pad_x = unit(1, "cm"), - pad_y = unit(1, "cm"), - height = unit(0.30, "cm"), - text_cex = 1 - ) + - annotation_north_arrow( - location = "br", - which_north = "true", - pad_x = unit(0.2, "cm"), - pad_y = unit(0.2, "cm") - ) + - theme( - axis.text = element_text(size = 12, family = "serif"), - plot.caption = element_text(size = 10, family = "serif"), - legend.text = element_text(size = 12, family = "serif"), - legend.title = element_text(size = 12, family = "serif"), - panel.background = element_rect(fill = "white"), - panel.grid.major = element_line(colour = "grey") - ) -``` - -Create 5km and 10km-resolution raster in NC shapefile - -```{r} -grid_10km <- rast(ext(nc_borders), resolution = c(10000, 10000), - crs = "EPSG:5070") -grid_10km <- extend(grid_10km, c(1, 1)) -grid_10km <- as.polygons(grid_10km) -grid_10km <- terra::intersect(grid_10km, nc_borders) -plot(grid_10km) - -grid_5km <- rast(ext(nc_borders), resolution = c(5000, 5000), crs = "EPSG:5070") -grid_5km <- extend(grid_5km, c(1, 1)) -grid_5km <- as.polygons(grid_5km) -grid_5km <- terra::intersect(grid_5km, nc_borders) -plot(grid_5km) -``` - -Compute and save 5kmx5km imperviousness file (each 5kmx5km cell is the mean of 30mx30m info) - -```{r, eval=F} -imp_mean <- zonal(imp, grid_5km, fun = "mean", as.raster = TRUE) -new_imp_file <- "NC_imperviousness_2019_5kmx5km.tif" -writeRaster(imp_mean, - filename = paste0(input_path, new_imp_file), - overwrite = TRUE) -``` - -Open and map the 5kmx5km imperviousness file - -```{r} -imp_mean <- rast("../input/NC_imperviousness_2019_5kmx5km.tif") -ggplot() + - geom_spatraster(data = imp_mean) + - geom_sf( - data = st_as_sf(nc_borders), aes(geometry = geometry), - colour = "white", linewidth = .3, fill = NA - ) + - scale_fill_whitebox_c( - palette = "muted", - labels = scales::label_number(suffix = "%"), - n.breaks = 12, - guide = guide_legend(reverse = TRUE) - ) + - labs( - fill = "", - title = "North Carolina imperviousness rate" - ) + - annotation_scale( - location = "bl", pad_x = unit(1, "cm"), - pad_y = unit(1, "cm"), - height = unit(0.30, "cm"), - text_cex = 1 - ) + - annotation_north_arrow( - location = "br", - which_north = "true", - pad_x = unit(0.2, "cm"), - pad_y = unit(0.2, "cm") - ) + - theme( - axis.text = element_text(size = 12, family = "serif"), - plot.caption = element_text(size = 10, family = "serif"), - legend.text = element_text(size = 12, family = "serif"), - legend.title = element_text(size = 12, family = "serif"), - panel.background = element_rect(fill = "white"), - panel.grid.major = element_line(colour = "grey") - ) -``` - -### Create urban and rural masks - -- Urban pixels: 5kmx5km imperviousness \> 5% - -- Urban pixels: 5kmx5km imperviousness \<= 5% - -```{r} -urb_mask <- ifel(imp_mean > 5, 1, NA) -rur_mask <- ifel(imp_mean <= 5, 1, NA) -ggplot() + - geom_spatraster(data = urb_mask) + - geom_sf( - data = st_as_sf(nc_borders), aes(geometry = geometry), - colour = "grey", linewidth = .3, fill = NA - ) + - scale_fill_whitebox_c( - palette = "deep", - labels = scales::label_number(suffix = ""), - na.value = NA - ) + - labs( - fill = "", - title = "North Carolina urban mask" - ) + - annotation_scale( - location = "bl", pad_x = unit(1, "cm"), - pad_y = unit(1, "cm"), - height = unit(0.30, "cm"), - text_cex = 1 - ) + - annotation_north_arrow( - location = "br", - which_north = "true", - pad_x = unit(0.2, "cm"), - pad_y = unit(0.2, "cm") - ) + - theme( - axis.text = element_text(size = 12, family = "serif"), - plot.caption = element_text(size = 10, family = "serif"), - legend.position = "none", - panel.background = element_rect(fill = "white"), - panel.grid.major = element_line(colour = "grey") - ) -``` - -Urban grid (200m-resolution) - -```{r} -urb_pol <- as.polygons(urb_mask) -urb_rast <- rast(urb_pol, resolution = c(200, 200), crs = "EPSG:5070") -urb_grid <- as.polygons(urb_rast) -urb_grid <- terra::intersect(urb_grid, urb_pol) -``` - -Rural grid (1km-resolution) - -```{r} -rur_pol <- as.polygons(rur_mask) -rur_rast <- rast(rur_pol, resolution = c(1000, 1000), crs = "EPSG:5070") -rur_grid <- as.polygons(rur_rast) -rur_grid <- terra::intersect(rur_grid, rur_pol) -``` - -### Map and save empty prediction grid - -```{r} -m <- ggplot() + - geom_spatvector(data = urb_grid, color = "orange", size = .1) + - geom_spatvector(data = rur_grid, color = "green", size = .1) + - geom_sf( - data = st_as_sf(nc_borders), aes(geometry = geometry), - colour = "grey", linewidth = .3, fill = NA - ) + - labs( - fill = "", - title = "Prediction grid" - ) + - annotation_scale( - location = "bl", pad_x = unit(1, "cm"), - pad_y = unit(1, "cm"), - height = unit(0.30, "cm"), - text_cex = 1 - ) + - annotation_north_arrow( - location = "br", - which_north = "true", - pad_x = unit(0.2, "cm"), - pad_y = unit(0.2, "cm") - ) + - theme( - axis.text = element_text(size = 12, family = "serif"), - plot.caption = element_text(size = 10, family = "serif"), - legend.text = element_text(size = 12, family = "serif"), - legend.title = element_text(size = 12, family = "serif"), - panel.background = element_rect(fill = "white"), - panel.grid.major = element_line(colour = "grey") - ) -m -``` - -Save urban and rural grid (squared polygons) and save all prediction grid as a unique dataframe (centroids of raster's polygons, coordinates in wgs84) - -```{r} -urb_grid_df <- centroids(urb_grid, inside = TRUE) -rur_grid_df <- centroids(rur_grid, inside = TRUE) -urb_grid_df$geo <- "urb" -rur_grid_df$geo <- "rur" -grid_points <- rbind(urb_grid_df, rur_grid_df) -grid_points <- project(grid_points, "EPSG:5070") -# -- keep only the geo covariate -grid_points <- grid_points %>% select(geo) -writeVector(grid_points, - paste0(input_path, "prediction_grid_points_urb_rur_empty.shp"), - overwrite = TRUE) -``` - -## Spatial covariates addition - -Interesting: with the terra::extract function, possibility to chose method='bilinear' to return a value interpolated from the values of the four nearest raster cells. - -```{r} -grid_points <- vect("../input/prediction_grid_points_urb_rur_empty.shp") -``` - -Create entent for RTP area if we want to plot a zoom - -```{r} -lat <- c(35.6, 36.11, 36.11, 35.6) -lon <- c(-79.19, -79.10, -78.39, -78.39) -ext_rtp <- vect(cbind(lon, lat), type = "points", crs = "EPSG:4326") -ext_rtp <- project(ext_rtp, crs(nc_borders)) -ext(ext_rtp) -``` - -#### Open covariates projected on "EPSG:5070" - -```{r} -imp_file <- "NC_imperviousness_2019_crs-meters.tif" -tcc_file <- "NC_tree-canopy-cover_2021_crs-meters.tif" -build_fp_file <- paste0("input/NC_building-footprints/", - "NorthCarolina_sum_crs-meters.tif") -build_h_file <- paste0("NC_building-height-by-block/", - "NC_building-heights-by-block_crs-meters.shp") -dem_file <- "NC-DEM_crs-meters.tif" - -imp <- rast(paste0(input_path, imp_file)) -tcc <- rast(paste0(input_path, tcc_file)) -build_fp <- rast(paste0(input_path, build_fp_file)) -build_h <- vect(paste0(input_path, build_h_file)) -dem <- rast(paste0(input_path, dem_file)) -``` - -#### Extract at each grid point - -```{r} -grid_spatcov <- grid_points -grid_spatcov <- terra::extract(imp, grid_points, bind = TRUE) -grid_spatcov <- terra::extract(tcc, grid_spatcov, - fun = function(x) mean(x, na.rm = TRUE), - method = "bilinear", bind = TRUE) -grid_spatcov <- terra::extract(build_fp, grid_spatcov, fun = "mean", - method = "bilinear", bind = TRUE) -grid_spatcov <- terra::extract(dem, grid_spatcov, bind = TRUE) - -# -- build_h is a vector and not a raster, bind=TRUE doesn't work... -grid0 <- terra::extract(build_h[, c("Height_cat")], grid_points) -grid_spatcov$build_h <- grid0$Height_cat -``` - -#### Rename covariates - -```{r} -grid_spatcov <- grid_spatcov %>% rename( - imp = "NC_imperviousness_2019", - tcc = "NC_tree-canopy-cover_2021", - build_fp = "NorthCarolina_sum", - dem = "Layer_1" -) - -summary(grid_spatcov) -``` - -## Temporal covariates addition - -TN with WMO standards - -```{r} -v <- data.frame(cbind(lon = geom(grid_points)[, "x"], - lat = geom(grid_points)[, "y"])) -grid_timecov <- vect(v, geom = c("lon", "lat")) -grid_timecov$geom <- geom(grid_timecov)[, "geom"] - -file_era5_tnwmo <- "era5_daily_reanalysis_20220601_20220831_TNwmo.tif" -file_era5_tn7am <- "era5_daily_reanalysis_20220601_20220831_TN7am.tif" -file_era5_tn12am <- "era5_daily_reanalysis_20220601_20220831_TN12am.tif" - -era5_tnwmo <- rast(paste0(input_path, file_era5_tnwmo)) -era5_tn7am <- rast(paste0(input_path, file_era5_tn7am)) -era5_tn12am <- rast(paste0(input_path, file_era5_tn12am)) - -era5_tnwmo <- project(era5_tnwmo, "EPSG:5070") -era5_tn7am <- project(era5_tn7am, "EPSG:5070") -era5_tn12am <- project(era5_tn12am, "EPSG:5070") - -grid_timecov_tnwmo <- terra::extract(era5_tnwmo, grid_timecov, bind = TRUE) -grid_timecov_tn7am <- terra::extract(era5_tn7am, grid_timecov, bind = TRUE) -grid_timecov_tn12am <- terra::extract(era5_tn12am, grid_timecov, bind = TRUE) - -file_era5_txwmo <- "era5_daily_reanalysis_20220601_20220831_TXwmo.tif" -file_era5_tx7am <- "era5_daily_reanalysis_20220601_20220831_TX7am.tif" -file_era5_tx12am <- "era5_daily_reanalysis_20220601_20220831_TX12am.tif" - -era5_txwmo <- rast(paste0(input_path, file_era5_txwmo)) -era5_tx7am <- rast(paste0(input_path, file_era5_tx7am)) -era5_tx12am <- rast(paste0(input_path, file_era5_tx12am)) - -era5_txwmo <- project(era5_txwmo, "EPSG:5070") -era5_tx7am <- project(era5_tx7am, "EPSG:5070") -era5_tx12am <- project(era5_tx12am, "EPSG:5070") - -grid_timecov_txwmo <- terra::extract(era5_txwmo, grid_timecov, bind = TRUE) -grid_timecov_tx7am <- terra::extract(era5_tx7am, grid_timecov, bind = TRUE) -grid_timecov_tx12am <- terra::extract(era5_tx12am, grid_timecov, bind = TRUE) - -df_grid_timecov_tnwmo <- melt(setDT(as.data.frame(grid_timecov_tnwmo)), - id.vars = "geom", - variable.name = "DATE", - value.name = "TNwmo") -df_grid_timecov_tn7am <- melt(setDT(as.data.frame(grid_timecov_tn7am)), - id.vars = "geom", - variable.name = "DATE", - value.name = "TN7am") -df_grid_timecov_tn12am <- melt(setDT(as.data.frame(grid_timecov_tn12am)), - id.vars = "geom", - variable.name = "DATE", - value.name = "TN12am") -df_grid_timecov_txwmo <- melt(setDT(as.data.frame(grid_timecov_txwmo)), - id.vars = "geom", - variable.name = "DATE", - value.name = "TXwmo") -df_grid_timecov_tx7am <- melt(setDT(as.data.frame(grid_timecov_tx7am)), - id.vars = "geom", - variable.name = "DATE", - value.name = "TX7am") -df_grid_timecov_tx12am <- melt(setDT(as.data.frame(grid_timecov_tx12am)), - id.vars = "geom", - variable.name = "DATE", - value.name = "TX12am") - -dt_time <- merge(df_grid_timecov_tnwmo, df_grid_timecov_tn7am, - by = c("geom", "DATE")) -dt_time <- merge(dt_time, df_grid_timecov_tn12am, by = c("geom", "DATE")) -dt_time <- merge(dt_time, df_grid_timecov_txwmo, by = c("geom", "DATE")) -dt_time <- merge(dt_time, df_grid_timecov_tx7am, by = c("geom", "DATE")) -dt_time <- merge(dt_time, df_grid_timecov_tx12am, by = c("geom", "DATE")) -``` - -## Merge spatial and temporal covariates - -Reproject grid points, add lat, lon, geom and transform to data.table - -```{r} -grid_spatcov <- project(grid_spatcov, "EPSG:4326") -grid_spatcov$geom <- geom(grid_spatcov)[, "geom"] -grid_spatcov$lon <- geom(grid_spatcov)[, "x"] -grid_spatcov$lat <- geom(grid_spatcov)[, "y"] -dt_spat <- setDT(as.data.frame(grid_spatcov)) -``` - -```{r} -grid_stcov <- merge(dt_spat, dt_time, by = c("geom")) -``` - -The prediction grid will be stored as a dataframe (or datatable) with columns: - -- lat, lon, date - -- TN, TX (varies in space and time) - -- dem, imp, tcc, build_fp, build_h (varies in space) - -It will have 56.794.820 rows (617335 pixels and 92 days) and 10 columns. - -```{r} -file <- "prediction_grid_points_urb_rur_space_time_covariates_jja2022.csv" -fwrite(grid_stcov, - paste0(input_path, file)) -``` - ------------------------------------------------------------------------- - -### Session info - -```{r} -sessionInfo() -``` diff --git a/vignettes/create_testdata.Rmd b/vignettes/create_testdata.Rmd index 94a5109..8260510 100755 --- a/vignettes/create_testdata.Rmd +++ b/vignettes/create_testdata.Rmd @@ -20,8 +20,8 @@ rtp_subset <- function(sp) { type = "polygons", crs = "EPSG:4326") ext_proj <- terra::project(ext, terra::crs(sp)) - crop_raster <- terra::crop(sp, ext_proj) - return(crop_raster) + crop_sp <- terra::crop(sp, ext_proj) + return(crop_sp) } ``` diff --git a/vignettes/divide-prediction-grid-per-day.Rmd b/vignettes/divide-prediction-grid-per-day.Rmd deleted file mode 100755 index 8bb16b4..0000000 --- a/vignettes/divide-prediction-grid-per-day.Rmd +++ /dev/null @@ -1,45 +0,0 @@ ---- -title: "divide-prediction-grid-per-day" -output: html_document -date: "2023-09-07" ---- - -```{r, message=F} -library(data.table) -``` - -Open full prediction grid - -```{r} -input_path <- "../input/" -pred <- fread(paste0(input_path, - "prediction_grid_points_", - "urb_rur_space_time_covariates_jja2022.csv")) -colnames(pred) <- c( - "geom", "geo", "imp", "tcc", "build.fp", - "dem", "build.h", "lon", "lat", "date", - "TNwmo", "TN7am", "TN12am", "TXwmo", "TX7am", "TX12am" -) -``` - -Create a folder to store all prediction grid per dates - -```{r} -if (!dir.exists("../input/prediction-grid")) { - dir.create("../input/prediction-grid") -} -``` - -Save each date seperately - -```{r} -dates <- seq(as.Date("2022-06-01"), as.Date("2022-08-31"), by = "1 day") -dates <- as.character(dates) -for (d in dates) { - file <- pred[date == d, ] - fwrite(file, paste0(input_path, - "prediction-grid/prediction_grid_points", - "_urb_rur_space_time_covariates_", d, ".csv")) -} -``` - diff --git a/vignettes/open-format-forest-height.Rmd b/vignettes/open-format-forest-height.Rmd deleted file mode 100755 index bf4e2cd..0000000 --- a/vignettes/open-format-forest-height.Rmd +++ /dev/null @@ -1,57 +0,0 @@ ---- -title: "open-format-forest-height" -output: html_document -date: "2023-08-30" ---- - -```{r, message=FALSE} -library(terra) -library(sf) -library(ggplot2) -library(ggspatial) -library(tidyterra) -library(rgeos) -library(data.table) -library(DT) -library(lubridate) -library(xts) -``` - -```{r} -input_path <- "../input/" -``` - -Open raw file of forest height on US - -```{r} -forest <- rast("../input/US_forest_height_2019.tif") -crs(forest) -``` - -Extract NC borders - -```{r} -nc_poly <- paste0("NC_county_boundary/", - "North_Carolina_State_and_County_Boundary_Polygons.shp") -nc_borders <- vect(paste0(input_path, nc_poly)) -nc_borders <- project(nc_borders, crs(forest)) -nc_forest <- terra::crop(forest, nc_borders) -``` - -Project on "EPSG:5070" - -```{r} -nc_forest_proj <- project(nc_forest, "EPSG:5070") -``` - -Save all final files - -```{r} -writeRaster(nc.forest, - filename = "../input/NC_forest_height_2019_crs-wgs84.tif", - overwrite = TRUE) -writeRaster(nc_forest_proj, - filename = paste0("../input/NC_forest_height_2019_albers.tif"), - overwrite = TRUE) -``` - diff --git a/vignettes/open-format-nlcd.Rmd b/vignettes/open-format-nlcd.Rmd deleted file mode 100755 index a96709d..0000000 --- a/vignettes/open-format-nlcd.Rmd +++ /dev/null @@ -1,65 +0,0 @@ ---- -title: "open-format-wb-height" -output: html_document -date: "2023-08-30" ---- - -```{r} -library(terra) -library(sf) -library(ggplot2) -library(ggspatial) -library(tidyterra) -library(rgeos) -library(data.table) -library(DT) -library(lubridate) -library(xts) -``` - - -Open raw file of land cover on US - -```{r} -input_path <- "../input/" -nlcd <- rast(paste0(input_path, - "/nlcd_2021_land_cover_l48_20230630/", - "nlcd_2021_land_cover_l48_20230630.img")) -``` - -Extract NC borders - -```{r} -nc_poly <- paste0("NC_county_boundary/", - "North_Carolina_State_and_County_Boundary_Polygons.shp") -nc_borders <- vect(paste0(input_path, nc_poly)) -nc_borders <- project(nc_borders, crs(nlcd)) -nc_nlcd <- terra::crop(nlcd, nc_borders) -``` - -Plot nc_nlcd - - -```{r} -plot(nc_nlcd) -plot(nc_borders, add = TRUE) -``` - - -Project on crs_proj - -```{r} -crs_proj <- "EPSG:5070" -nc_nlcd_proj <- project(nc_nlcd, crs_proj) -``` - -Save all final files - -```{r} -writeRaster(nc_nlcd, - filename = paste0(input_path, "NC_nlcd_crs-wgs84.tif"), - overwrite = TRUE) -writeRaster(nc_nlcd_proj, - filename = paste0("../input/NC_nlcd_crs-albers.tif"), - overwrite = TRUE) -``` diff --git a/vignettes/open-project-spatial-covariates.Rmd b/vignettes/open-project-spatial-covariates.Rmd deleted file mode 100755 index ea27c55..0000000 --- a/vignettes/open-project-spatial-covariates.Rmd +++ /dev/null @@ -1,122 +0,0 @@ ---- -title: "open-project-spatial-covariates.R" -output: html_document -date: "2023-08-22" ---- - -Libraries - -```{r, message=F} -library(terra) -library(sf) -library(ggplot2) -library(ggspatial) -library(tidyterra) -library(rgeos) -library(data.table) # -- for large flat datasets -``` - -```{r} -input_path <- "../input/" -``` - -Open data - -```{r} -imp_file <- "NC_imperviousness_2019_crs-meters.tif" -tcc_file <- "NC_tree-canopy-cover_2021_crs-meters.tif" -build_fp_file <- paste0("input/NC_building-footprints/", - "NorthCarolina_sum_crs-meters.tif") -build_h_file <- paste0("NC_building-height-by-block/", - "NC_building-heights-by-block_crs-meters.shp") -dem_file <- "NC-DEM_crs-meters.tif" - -imp <- rast(paste0(input_path, imp_file)) -tcc <- rast(paste0(input_path, tcc_file)) -build_fp <- rast(paste0(input_path, build_fp_file)) -build_h <- vect(paste0(input_path, build_h_file)) -dem <- rast(paste0(input_path, dem_file)) -``` - -```{r} -crs(imp) -cat("\n") -crs(tcc) -cat("\n") -crs(build_fp) -cat("\n") -crs(build_h) -cat("\n") -crs(dem) -``` - -Choose a projection - -```{r} -crs <- "EPSG:4326" -crs_name <- "crs-wgs84" -cat(crs) -``` - -Project data and save the data - -```{r} -# -- imperviousness -imp_p <- project(imp, crs) -writeRaster(imp_p, - filename = paste0(input_path, - "NC_imperviousness_2019_", - crs_name, - ".tif"), - overwrite = TRUE) - -# -- tree canopy cover -tcc_p <- project(tcc, crs) -writeRaster(tcc_p, - filename = paste0(input_path, - "NC_tree-canopy-cover_2021_", - crs_name, - ".tif"), - overwrite = TRUE) - -# -- building footprint -build_fp_p <- project(build_fp, crs) -writeRaster(build_fp_p, - filename = paste0(input_path, - "NC_building-footprints/NorthCarolina_sum_", - crs_name, - ".tif"), - overwrite = TRUE) - -# -- building height -build_h_p <- project(build_h, crs) -writeVector(build_h_p, - filename = paste0(input_path, - "NC_building-height-by-block/", - "NC_building-heights-by-block_", - crs_name, - ".shp"), - overwrite = TRUE) - -# -- digital elevation model -dem_p <- project(dem, crs) -writeRaster(dem_p, - filename = paste0(input_path, - "NC-DEM_", - crs_name, - ".tif"), - overwrite = TRUE) -``` - -Important: it sounds like there are missing tiles in building footprint data. -An explanation can be found in the paper Heris, M.P., Foks, N., Bagstad, K., -and Troy, A., 2020, A national dataset of rasterized building footprints for -the U.S.: U.S. Geological Survey data release, -: - -"*We also identified systematic gaps in the Microsoft data for some geographic -areas. These larger gaps seem to have a tile pattern, where aerial photos may -have been unavailable to the Microsoft building detection algorithm"* - -Their computational algorithm is applied to Microsoft released a U.S.-wide -vector building dataset provided in 2018 but this dataset has missing tiles.