Skip to content

Commit

Permalink
add testdata loading file function and fix create_prediction_grid pro…
Browse files Browse the repository at this point in the history
…blems (especially add_era5_raster)
  • Loading branch information
Marques committed Dec 18, 2023
1 parent 641e159 commit dfcbbfa
Show file tree
Hide file tree
Showing 9 changed files with 93 additions and 87 deletions.
49 changes: 25 additions & 24 deletions R/add_covariates.R
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
#' @returns the same terra::SpatVector with the column "imp"
#' @importFrom magrittr "%>%"
#' @export
add_imp <- function(imp_path, sp_vect) {
add_imp <- function(sp_vect, imp_path) {
if (!file.exists(imp_path)) {
stop("imp_path does not exist")
}
Expand All @@ -15,7 +15,7 @@ add_imp <- function(imp_path, sp_vect) {
dplyr::rename(imp = names(imp)) %>%
terra::project(terra::crs(sp_vect))
if (any(is.na(sp_vect_cov$imp))) {
stop("NAs found in imp column")
warning("NAs found in imp column")
}
return(sp_vect_cov)
}
Expand All @@ -27,7 +27,7 @@ add_imp <- function(imp_path, sp_vect) {
#' @param sp_vect a terra::SpatVector
#' @returns the same terra::SpatVector with the column "tcc"
#' @export
add_tcc <- function(tcc_path, sp_vect) {
add_tcc <- function(sp_vect, tcc_path) {
if (!file.exists(tcc_path)) {
stop("tcc_path does not exist")
}
Expand All @@ -41,7 +41,7 @@ add_tcc <- function(tcc_path, sp_vect) {
dplyr::rename(tcc = names(tcc)) %>%
terra::project(terra::crs(sp_vect))
if (any(is.na(sp_vect_cov$tcc))) {
stop("NAs found in tcc column")
warning("NAs found in tcc column")
}
return(sp_vect_cov)
}
Expand All @@ -52,7 +52,7 @@ add_tcc <- function(tcc_path, sp_vect) {
#' @param sp_vect a terra::SpatVector
#' @returns the same terra::SpatVector with the column "dem"
#' @export
add_dem <- function(dem_path, sp_vect) {
add_dem <- function(sp_vect, dem_path) {
if (!file.exists(dem_path)) {
stop("dem_path does not exist")
}
Expand All @@ -62,7 +62,7 @@ add_dem <- function(dem_path, sp_vect) {
dplyr::rename(dem = names(dem)) %>%
terra::project(terra::crs(sp_vect))
if (any(is.na(sp_vect_cov$dem))) {
stop("NAs found in dem column")
warning("NAs found in dem column")
}
return(sp_vect_cov)
}
Expand All @@ -74,7 +74,7 @@ add_dem <- function(dem_path, sp_vect) {
#' @returns the same terra::SpatVector with terrain covariates
#' (dem, slope, aspect, roughness, flowdir)
#' @export
add_terrain <- function(dem_path, sp_vect) {
add_terrain <- function(sp_vect, dem_path) {
dem_rast <- terra::rast(dem_path)
names(dem_rast) <- "dem"
dem_rast$slope <- terra::terrain(dem_rast$dem, "slope")
Expand All @@ -95,19 +95,19 @@ add_terrain <- function(dem_path, sp_vect) {
terra::extract(x = dem_rast, bind = TRUE) %>%
terra::project(terra::crs(sp_vect))
if (any(is.na(sp_vect_cov$dem))) {
stop("NAs found in dem column")
warning("NAs found in dem column")
}
if (any(is.na(sp_vect_cov$slope))) {
stop("NAs found in slope column")
warning("NAs found in slope column")
}
if (any(is.na(sp_vect_cov$aspect))) {
stop("NAs found in aspect column")
warning("NAs found in aspect column")
}
if (any(is.na(sp_vect_cov$roughness))) {
stop("NAs found in roughness column")
warning("NAs found in roughness column")
}
if (any(is.na(sp_vect_cov$flowdir))) {
stop("NAs found in flowdir column")
warning("NAs found in flowdir column")
}
return(sp_vect_cov)
}
Expand All @@ -121,7 +121,7 @@ add_terrain <- function(dem_path, sp_vect) {
#' @param sp_vect a terra::SpatVector
#' @returns the same terra::SpatVector with the column "build_fp"
#' @export
add_canopy_h <- function(canopy_h_path, sp_vect) {
add_canopy_h <- function(sp_vect, canopy_h_path) {
if (!file.exists(canopy_h_path)) {
stop("canopy_h_path does not exist")
}
Expand All @@ -137,7 +137,7 @@ add_canopy_h <- function(canopy_h_path, sp_vect) {
dplyr::rename(canopy_h = names(canopy_h)) %>%
terra::project(terra::crs(sp_vect))
if (any(is.na(sp_vect_cov$canopy_h))) {
stop("NAs found in canopy_h column")
warning("NAs found in canopy_h column")
}
return(sp_vect_cov)
}
Expand All @@ -150,7 +150,7 @@ add_canopy_h <- function(canopy_h_path, sp_vect) {
#' @param sp_vect a terra::SpatVector
#' @returns the same terra::SpatVector with the column "build_fp"
#' @export
add_build_fp <- function(build_fp_path, sp_vect) {
add_build_fp <- function(sp_vect, build_fp_path) {
if (!file.exists(build_fp_path)) {
stop("build_fp_path does not exist")
}
Expand All @@ -165,7 +165,7 @@ add_build_fp <- function(build_fp_path, sp_vect) {
dplyr::rename(build_fp = names(build_fp)) %>%
terra::project(terra::crs(sp_vect))
if (any(is.na(sp_vect_cov$build_fp))) {
stop("NAs found in build_fp column")
warning("NAs found in build_fp column")
}
return(sp_vect_cov)
}
Expand All @@ -177,8 +177,7 @@ add_build_fp <- function(build_fp_path, sp_vect) {
#' @param sp_vect a terra::SpatVector
#' @returns the same terra::SpatVector with the column "build_h"
#' @export
add_build_h <- function(
build_h_path, sp_vect) {
add_build_h <- function(sp_vect, build_h_path) {
if (!file.exists(build_h_path)) {
stop("build_h_path does not exist")
}
Expand All @@ -189,7 +188,7 @@ add_build_h <- function(
sp_vect_cov$build_h <- sp_vect_build_h$Height_cat
sp_vect_cov <- terra::project(sp_vect_cov, terra::crs(sp_vect))
if (all(is.na(sp_vect_cov$build_h) | all(sp_vect_cov$build_h == "NA"))) {
stop("build_h column is only NA")
warning("build_h column is only NA")
}
return(sp_vect_cov)
}
Expand Down Expand Up @@ -267,7 +266,7 @@ add_nlcd_ratio <- function(data_vect,
#' @param datatable A data.table object with columns "lat", "lon"
#' @param crs A character containing the crs of spatial data
#' @returns same datatable object with "county" columns
add_county <- function(county_path, data_vect) {
add_county <- function(data_vect, county_path) {
if (!file.exists(county_path)) {
stop("county_path does not exist")
}
Expand Down Expand Up @@ -373,29 +372,31 @@ add_era5_vect <- function(data_vect, era5_path) {
#' @param era5_path A character path to era5 csv
#' @returns A SpatRasterDataset
#' @export
add_era5_rast <- function(data_rast, era5_path = era5) {
add_era5_rast <- function(data_rast, era5_path) {
era5 <- data.table::fread(era5_path) %>%
HeatModel::create_stdtobj(crs_stdt = "EPSG:4326") %>%
HeatModel::convert_stdt_spatrastdataset()
data_vect <- terra::as.points(data_rast)
# empty prediction SpatVector
new_data_vect <- terra::vect(terra::geom(data_vect)[, c("x", "y")],
new_data_vect <- terra::vect(
terra::geom(data_vect)[, c("x", "y")],
type = "points",
crs = terra::crs(data_vect)
)
# extract each daily covariate based on era5 and convert to raster
pred_rds_era5 <- list()
for (i in 2:7) {
names <- names(era5[[i]])
pred_rds_era5[[i - 1]] <- terra::project(
new_data_vect,
terra::crs(era5[[i]])
) %>%
terra::extract(x = era5[[i]], bind = TRUE) %>%
terra::project(terra::crs(new_data_vect)) %>%
terra::rasterize(data_rast, field = names())
# maybe names() won't work
terra::rasterize(data_rast, field = names)
}
# turn into a SpatRasterDataset
pred_rds_era5 <- terra::sds(pred_rds_era5)
names(pred_rds_era5) <- names(era5[[2:7]])
return(pred_rds_era5)
}
63 changes: 18 additions & 45 deletions R/create_prediction_grid_300m_covariates.R
Original file line number Diff line number Diff line change
@@ -1,44 +1,40 @@
#' Create prediction grid with the help of imperviousness raster file
#'
#' @param borders_path a character path to borders shapefile
#' @param covar_path a list of character path to covariates files
#' @param covar_files_list a list of character path to covariates files
#' @param output_path a character path to the output storage folder
#' @returns the same terra::SpatVector with the column "imp"
#' @export
create_pred_rds <- function(borders_path,
covar_path,
covar_files_list,
output_path) {
# create spat raster in borders
borders <- terra::vect(borders_path)
imp <- terra::rast(covar_path$imp)
imp <- terra::rast(covar_files_list$imp)
borders_proj <- terra::project(borders, terra::crs(imp))
imp <- terra::mask(borders_proj)
imp <- terra::mask(imp, borders_proj)

# aggregation at 300m
pred_rast <- aggregate(imp, fact = 10, fun = "mean")
pred_rast <- terra::aggregate(imp, fact = 10, fun = "mean")
names(pred_rast) <- "imp"

# vector format will also be useful for covariate extraction
pred_vect <- terra::as.points(pred_rast)

# ADD SPATIAL COVARIATES
pred_vect <- add_canopy_h(covar_path$canopy_h, sp_vect = pred_vect) %>%
add_terrain(covar_path$dem) %>%
add_tcc(covar_path$dem) %>%
add_build_fp(covar_path$build_fp) %>%
add_build_h(covar_path$build_h) %>%
add_canopy_h(covar_path$canopy_h) %>%
pred_vect <- add_canopy_h(covar_files_list$canopy_h, sp_vect = pred_vect) %>%
add_terrain(covar_files_list$dem) %>%
add_tcc(covar_files_list$tcc) %>%
add_build_fp(covar_files_list$build_fp) %>%
add_build_h(covar_files_list$build_h) %>%
add_nlcd_ratio(buf_radius = 150,
nlcd_path = covar_path$nlcd) %>%
add_county(county_path = covar_path$county)

nlcd_path = covar_files_list$nlcd) %>%
add_county(county_path = covar_files_list$county)
# turn into a raster
pred_rast <- terra::rasterize(pred_vect,
pred_rast,
field = names(pred_vect)
)
# Terrain covariates
pred_rast <- add_terrain(pred_rast)
# Create SpatRasterDataset with all spatial covariates
pred_rds <- list()
for (i in names(pred_rast)) {
Expand All @@ -47,45 +43,22 @@ create_pred_rds <- function(borders_path,
pred_rds <- terra::sds(pred_rds)

# ADD ERA5 COVARIATES
era5 <- fread(covar_path$era5) %>%
data.table::rename(time = date)
# convert era5 to stdt and then to SpatRasterDataset
era5_stdt <- create_stdtobj(era5, "EPSG:4326")
era5_rds <- convert_stdt_spatrastdataset(era5_stdt)
# empty prediction SpatVector
new_pred_vect <- terra::vect(terra::geom(pred_vect)[, c("x", "y")],
type = "points",
crs = terra::crs(pred_vect)
)
# extract each daily covariate based on era5 and convert to raster
pred_rds_era5 <- list()
for (i in 2:7) {
pred_rds_era5[[i - 1]] <- terra::project(
new_pred_vect,
terra::crs(era5_rds[[i]])
) %>%
terra::extract(x = era5_rds[[i]], bind = TRUE) %>%
terra::project(terra::crs(new_pred_vect)) %>%
terra::rasterize(pred_rast, field = terra::names())
}
# turn into a SpatRasterDataset
pred_rds_era5 <- terra::sds(pred_rds_era5)
terra::names(pred_rds_era5) <- terra::names(era5_rds[[2:7]])

pred_rds_era5 <- add_era5_rast(pred_rast, era5_path = covar_files_list$era5)
# SAVE OUTPUT PER DAY AS A .TIF
dates <- terra::names(pred_rds_era5[[1]])
era5_covs <- terra::names(pred_rds_era5)
dates <- names(pred_rds_era5[[1]])
era5_covs <- names(pred_rds_era5)
list_rast_dates <- list()
for (date in dates) {
rast_date <- pred_rast
for (era5_cov in era5_covs) {
rast_date[[era5_cov]] <- pred_rds_era5[era5_cov][date]
rast_date[[era5_cov]] <- pred_rds_era5[[era5_cov]][[date]]
}
list_rast_dates[[date]] <- rast_date
terra::writeRaster(rast_date,
filename = paste0(output_path,
"pred_grid_",
date,
".tif"))
".tif"),
overwrite = TRUE)
}
}
23 changes: 23 additions & 0 deletions R/upload_metadata.R
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,29 @@ list_covar_nc <- function(covar_folder) {
}


#' Give a list of covariate files names
#'
#' @param covar_folder character path to covariates files
#' @returns a list of covariate file paths
#' @export
list_covar_testdata <- function(covar_folder) {
covar_files <- list(
imp = paste0(covar_folder, "rtp_imp.tif"),
dem = paste0(covar_folder, "rtp_dem.tif"),
tcc = paste0(covar_folder, "rtp_tcc.tif"),
canopy_h = paste0(covar_folder, "rtp_canopy_h.tif"),
build_h = paste0(covar_folder, "rtp_build_h.shp"),
build_fp = paste0(covar_folder, "rtp_build_fp.tif"),
nlcd = paste0(covar_folder, "rtp_nlcd.tif"),
county = paste0(covar_folder, "rtp_counties.shp"),
era5 = paste0(covar_folder, "rtp_era5.csv")
)
if (!all(sapply(covar_files, file.exists))) {
warning("Some of the files do not exist.")
}
return(covar_files)
}

#' Give a list of monitor files names
#'
#' @param folder character path to monitors folder
Expand Down
Binary file modified tests/testdata/rtp_spatvector.dbf
Binary file not shown.
Binary file modified tests/testdata/rtp_spatvector.shp
Binary file not shown.
Binary file modified tests/testdata/rtp_spatvector.shx
Binary file not shown.
Loading

0 comments on commit dfcbbfa

Please sign in to comment.