Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Several bug fixes #124

Merged
merged 24 commits into from
Dec 13, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions Example_Use.R
Original file line number Diff line number Diff line change
Expand Up @@ -84,11 +84,11 @@ ggplot(res2[RUN != "ensembleMean" & GCM != "Observed",],aes(x = PERIOD, y = PPT,
thebb <- get_bb(in_xyz) ##get bounding box based on input points
dbCon <- data_connect()
##get normal
normalbc <- normal_input_postgis(dbCon = dbCon, normal = "normal_bc", bbox = thebb, cache = TRUE)
normalbc <- normal_input(dbCon = dbCon, normal = "normal_bc", bbox = thebb, cache = TRUE)
plot(normalbc[[13]])

##gcm annomalies
gcm <- gcm_input_postgis(dbCon, bbox = thebb, gcm = c("ACCESS-ESM1-5", "EC-Earth3"),
gcm <- gcm_input(dbCon, bbox = thebb, gcm = c("ACCESS-ESM1-5", "EC-Earth3"),
ssp = c("ssp370"),
period = c("2021_2040","2041_2060","2061_2080"),
max_run = 0,
Expand Down
4 changes: 2 additions & 2 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ export(data_prepare)
export(data_update)
export(downscale)
export(gcm_hist_input)
export(gcm_input_postgis)
export(gcm_input)
export(gcm_ts_input)
export(get_bb)
export(historic_input)
Expand All @@ -24,7 +24,6 @@ export(list_run)
export(list_ssp)
export(list_variables)
export(normal_input)
export(normal_input_postgis)
export(pgGetTerra)
import(data.table)
import(uuid)
Expand Down Expand Up @@ -69,6 +68,7 @@ importFrom(terra,nlyr)
importFrom(terra,rast)
importFrom(terra,resample)
importFrom(terra,sources)
importFrom(terra,sprc)
importFrom(terra,writeRaster)
importFrom(terra,xres)
importFrom(terra,yres)
Expand Down
1 change: 1 addition & 0 deletions R/data.R
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
#' Update external package data
#' @param gcm A character. Relative path from the source root to global circulation models files folder.
#' Default to option value "climr.gcm.path" if set, or "inputs_pkg/gcmData".
#' @param historic blah blah
#' @param normal A character. Relative path from the source root to base normal files folder.
#' @param quiet A logical. If `TRUE`, suppress status messages (if any), and the progress bar.
#' @param ... Others parameters such as `source` or `repo` for content getting functions.
Expand Down
59 changes: 32 additions & 27 deletions R/dbGetRaster.R
Original file line number Diff line number Diff line change
@@ -1,34 +1,31 @@
#' Download raster with bounding box from PostGIS
#' @param conn a DBI or RPostgres connection object
#' @template conn
#' @param name character. Table name in database
#' @param tile TODO
#' @param rast character. Name of column which stores raster data.
#' Defaults to "rast"
#' @param bands Which raster bands to return. Default 37:73.
#' @template bands
#' @template boundary

#' @return A SpatRaster
#'
#' @importFrom data.table setDT copy
#' @importFrom terra rast merge
#' @importFrom terra rast merge sprc
#' @importFrom RPostgres dbQuoteIdentifier dbGetQuery
#' @importFrom DBI dbQuoteIdentifier dbGetQuery
#' @export

pgGetTerra <- function(conn, name, tile, rast = "rast", bands = 37:73,
boundary = NULL) {
## Check and prepare the schema.name
name1 <- name
nameque <- paste(name1, collapse = ".")
namechar <- gsub("'", "''", paste(gsub('^"|"$', "", name1), collapse = "."))

## rast query name
rastque <- dbQuoteIdentifier(conn, rast)

projID <- dbGetQuery(conn, paste0("select ST_SRID(", rastque, ") as srid from ", nameque, " where rid = 1;"))$srid[1]




if (length(bands) > 1664) { ## maximum number of columns
info <- dbGetQuery(conn, paste0(
"select
Expand All @@ -52,7 +49,7 @@ pgGetTerra <- function(conn, name, tile, rast = "rast", bands = 37:73,
bands_temp <- bands[brks[i]:(brks[i + 1] - 1)]
bandqs1 <- paste0("UNNEST(ST_Dumpvalues(rast, ", bands_temp, ")) as vals_", bands_temp)
bandqs2 <- paste0("ST_Union(rast", rastque, ",", bands_temp, ") rast_", bands_temp)

rast_vals_temp <- dbGetQuery(conn, paste0(
"SELECT ", paste(bandqs1, collapse = ","),
" from (SELECT ST_Union(rast) rast FROM ", nameque, " WHERE ST_Intersects(",
Expand All @@ -77,14 +74,14 @@ pgGetTerra <- function(conn, name, tile, rast = "rast", bands = 37:73,
)
} else {
if (!tile) {
rout <- make_raster(boundary)
rout <- make_raster(boundary, conn, rastque, nameque, projID, bands)
return(rout)
}
max_dist <- 5
# if(boundary[1] - boundary[2] > max_dist | boundary[3] - boundary[4] > max_dist) {
x_seq <- unique(c(seq(boundary[2], boundary[1], by = max_dist), boundary[1]))
y_seq <- unique(c(seq(boundary[4], boundary[3], by = max_dist), boundary[3]))

boundary_ls <- list()
if (length(x_seq) < 2 | length(y_seq) < 2) {
boundary_ls[["11"]] <- boundary
Expand All @@ -95,9 +92,12 @@ pgGetTerra <- function(conn, name, tile, rast = "rast", bands = 37:73,
}
}
}


r_list <- lapply(boundary_ls, FUN = make_raster)


r_list <- lapply(boundary_ls, FUN = make_raster,
conn = conn, rastque = rastque,
nameque = nameque, projID = projID,
bands = bands)
r_list <- r_list[!sapply(r_list, is.null)]
if (length(r_list) > 1) {
rout <- merge(sprc(r_list))
Expand Down Expand Up @@ -178,31 +178,36 @@ get_bb <- function(in_xyz) {
}



#' Make raster from a boundary
#'
#' Used internally to access the PostGRS database and
#' create a SpatRaster using a given spatial boundary
#'
#' @template boundary
#' @template boundary
#' @template conn
#' @param rastque rast query name obtained with e.g. `dbQuoteIdentifier(conn, "rast")`
#' @param nameque schema.name
#' @param projID projID in data.base
#' @template bands
#'
#' @return a SpatRaster
#'
#' @importFrom DBI dbGetQuery
#' @importFrom terra rast
#'
make_raster <- function(boundary) {
make_raster <- function(boundary, conn, rastque, nameque, projID, bands) {
cat(".")
info <- dbGetQuery(conn, paste0(
"select
st_xmax(st_envelope(rast)) as xmx,
st_xmin(st_envelope(rast)) as xmn,
st_ymax(st_envelope(rast)) as ymx,
st_ymin(st_envelope(rast)) as ymn,
st_width(rast) as cols,
st_height(rast) as rows
from
(select st_union(", rastque, ",", 1, ") rast from ", nameque, "\n
WHERE ST_Intersects(",
st_xmax(st_envelope(rast)) as xmx,
st_xmin(st_envelope(rast)) as xmn,
st_ymax(st_envelope(rast)) as ymx,
st_ymin(st_envelope(rast)) as ymn,
st_width(rast) as cols,
st_height(rast) as rows
from
(select st_union(", rastque, ",", 1, ") rast from ", nameque, "\n
WHERE ST_Intersects(",
rastque, ",ST_SetSRID(ST_GeomFromText('POLYGON((", boundary[4],
" ", boundary[1], ",", boundary[4], " ", boundary[2],
",\n ", boundary[3], " ", boundary[2], ",", boundary[3],
Expand Down
45 changes: 33 additions & 12 deletions R/downscale.R
Original file line number Diff line number Diff line change
@@ -1,12 +1,15 @@
#' Downscale and Calculate climate variables for points of interest

#'
#' @details
#' Additional details... TODO.
#'
#' @param xyz three or four column data.frame: long, lat, elev, (id)
#' @param which_normal Select which normal layer to use. Default is "auto", which selects the highest resolution normal for each point
#' @param historic_period Which historic period. Default `NULL`
#' @param historic_ts Historic years requested. Must be in `1902:2015`. Default `NULL`
#' @param gcm_models Vector of GCM names. Options are `list_gcm()`. Used for gcm periods, gcm timeseries, and historic timeseries. Default `NULL`
#' @template ssp
#' @param gcm_period Requested future periods. Options are `list_period()`
#' @param gcm_period Requested future periods. Options are `list_gcm_period()`
#' @param gcm_ts_years Requested future timeseries years. Must be in `2015:2100`
#' @param gcm_hist_years Requested historic modelled years. Must be in `1851:2010`
#' @template max_run
Expand Down Expand Up @@ -58,9 +61,9 @@ climr_downscale <- function(xyz, which_normal = c("auto", "BC", "NorAm"), histor

message("Getting normals...")
if (which_normal == "NorAm") {
normal <- normal_input_postgis(dbCon = dbCon, normal = "normal_na", bbox = thebb, cache = cache)
normal <- normal_input(dbCon = dbCon, normal = "normal_na", bbox = thebb, cache = cache)
} else if (which_normal == "BC") {
normal <- normal_input_postgis(dbCon = dbCon, normal = "normal_bc", bbox = thebb, cache = cache)
normal <- normal_input(dbCon = dbCon, normal = "normal_bc", bbox = thebb, cache = cache)
} else {
if (ncol(xyz) == 3) xyz$ID <- 1:nrow(xyz)
bc_outline <- rast(system.file("extdata", "bc_outline.tif", package = "climr"))
Expand All @@ -71,9 +74,9 @@ climr_downscale <- function(xyz, which_normal = c("auto", "BC", "NorAm"), histor
xyz <- xyz[!is.na(pnts$PPT01),]
thebb_bc <- get_bb(xyz)
message("for BC...")
normal <- normal_input_postgis(dbCon = dbCon, normal = "normal_bc", bbox = thebb_bc, cache = cache)
normal <- normal_input(dbCon = dbCon, normal = "normal_bc", bbox = thebb_bc, cache = cache)
} else {
normal <- normal_input_postgis(dbCon = dbCon, normal = "normal_na", bbox = thebb, cache = cache)
normal <- normal_input(dbCon = dbCon, normal = "normal_na", bbox = thebb, cache = cache)
}
}

Expand All @@ -87,7 +90,7 @@ climr_downscale <- function(xyz, which_normal = c("auto", "BC", "NorAm"), histor
if (!is.null(gcm_models)) {
if (!is.null(gcm_period)) {
message("Getting GCMs...")
gcm <- gcm_input_postgis(dbCon,
gcm <- gcm_input(dbCon,
bbox = thebb, gcm = gcm_models,
ssp = ssp,
period = gcm_period,
Expand Down Expand Up @@ -145,7 +148,7 @@ climr_downscale <- function(xyz, which_normal = c("auto", "BC", "NorAm"), histor
na_xyz <- xyz_save[!xyz_save[, 4] %in% bc_ids,]
thebb <- get_bb(na_xyz)
message("Now North America...")
normal <- normal_input_postgis(dbCon = dbCon, normal = "normal_na", bbox = thebb, cache = cache)
normal <- normal_input(dbCon = dbCon, normal = "normal_na", bbox = thebb, cache = cache)

results_na <- downscale(
xyz = na_xyz,
Expand All @@ -169,6 +172,10 @@ climr_downscale <- function(xyz, which_normal = c("auto", "BC", "NorAm"), histor


#' Downscale target rasters to points of interest
#'
#' @details
#' Additional details... TODO.
#'
#' @param xyz A 3-column matrix or data.frame (x, y, z) or (lon, lat, elev).
#' @param normal Reference normal baseline input from `normal_input`.
#' @param gcm Global Circulation Models input from `gcm_input`. Default to NULL.
Expand All @@ -182,21 +189,35 @@ climr_downscale <- function(xyz, which_normal = c("auto", "BC", "NorAm"), histor
#' `variables` dataset. Default to monthly PPT, Tmax, Tmin.
#' @param ppt_lr A boolean. Apply lapse rate adjustment to precipitations. Default to FALSE.
#' @param nthread An integer. Number of parallel threads to use to do computations. Default to 1L.
#'
#' @import data.table
#' @importFrom terra extract rast sources ext xres yres crop
#' @importFrom parallel makeForkCluster makePSOCKcluster stopCluster splitIndices parLapply
#'
#' @return A downscaled dataset. If `gcm` is NULL, this is just the downscaled `normal`
#' at point locations. If `gcm` is provided, this returns a downscaled dataset for each
#' point location, general circulation model, shared socioeconomic pathway, run and period.
#'
#' @export
#' @examples
#' \dontrun{
#' dbCon <- data_connect()
#' xyz <- data.frame(lon = runif(10, -140, -106), lat = runif(10, 37, 61), elev = runif(10))
#' normal <- normal_input()
#' gcm_input <- gcm_input(list_gcm()[3], list_ssp()[1], list_period()[2])
#'
#' ## get bounding box based on input points
#' thebb <- get_bb(xyz)
#' normal <- normal_input(dbCon = dbCon, bbox = thebb, cache = TRUE)
#'
#' ## pick one GCM, one SSP and one period from the list of available options
#' gcm <- gcm_input(dbCon, thebb, gcm = list_gcm()[3], list_ssp()[1], list_gcm_period()[2])
#'
#' ## notice coarseness of the data
#' terra::plot(gcm[[1]])
#'
#' downscale(xyz, normal, gcm)
#' historic <- historic_input()
#' out <- downscale(xyz, normal, gcm = NULL, historic = historic, ppt_lr = FALSE)
#' historic <- historic_input(dbCon, thebb)
#'
#' downscale(xyz, normal, gcm = NULL, historic = historic, ppt_lr = FALSE)
#' }
downscale <- function(xyz, normal, gcm = NULL, historic = NULL, gcm_ts = NULL, gcm_hist = NULL, historic_ts = NULL, return_normal = FALSE,
vars = sort(sprintf(c("PPT%02d", "Tmax%02d", "Tmin%02d"), sort(rep(1:12, 3)))),
Expand Down
11 changes: 7 additions & 4 deletions R/gcm.R
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@
#' @importFrom uuid UUIDgenerate
#' @import data.table
#' @export
gcm_input_postgis <- function(dbCon, bbox = NULL, gcm = list_gcm(), ssp = list_ssp(), period = list_period(), max_run = 0L, cache = TRUE) {
gcm_input <- function(dbCon, bbox = NULL, gcm = list_gcm(), ssp = list_ssp(), period = list_gcm_period(), max_run = 0L, cache = TRUE) {
dbnames <- structure(list(
GCM = c(
"ACCESS-ESM1-5", "BCC-CSM2-MR", "CanESM5",
Expand Down Expand Up @@ -107,7 +107,7 @@ gcm_input_postgis <- function(dbCon, bbox = NULL, gcm = list_gcm(), ssp = list_s
#' @template bbox
#' @template gcm
#' @param years Numeric vector of desired years. Must be in `1851:2015`.
#' Can be obtained from `list_period()`. Default to `list_period()`.
#' Can be obtained from `list_gcm_period()`. Default to `list_gcm_period()`.
#' @template max_run
#' @param cache Logical specifying whether to cache new data locally or no. Default `TRUE`
#' @return An object to use with `downscale`. A list of `SpatRaster` with, possibly, multiple layers.
Expand Down Expand Up @@ -205,21 +205,24 @@ gcm_hist_input <- function(dbCon, bbox = NULL, gcm = list_gcm(), years = 1901:19
# period <- 2020:2050

#' Create gcm timeseries input for `downscale` using data on Postgis database.
#'
#' @template dbCon
#' @template bbox
#' @template gcm
#' @template ssp
#' @param years Numeric or character vector in in `2020:2100`
#' @param years Numeric or character vector in `2020:2100`. Defaults to `2020:2030`.
#' @template max_run
#' @param cache Logical specifying whether to cache new data locally or no. Default `TRUE`
#'
#' @return An object to use with `downscale`. A list of `SpatRaster` with, possibly, multiple layers.
#'
#' @importFrom terra rast writeRaster ext nlyr
#' @importFrom utils head
#' @importFrom RPostgres dbGetQuery
#' @import uuid
#' @import data.table
#' @export
gcm_ts_input <- function(dbCon, bbox = NULL, gcm = list_gcm_ts(), ssp = list_ssp(), years = list_years_ts(), max_run = 0L, cache = TRUE) {
gcm_ts_input <- function(dbCon, bbox = NULL, gcm = list_gcm(), ssp = list_ssp(), years = 2020:2030, max_run = 0L, cache = TRUE) {
period <- years
dbnames <- structure(list(
GCM = c(
Expand Down
6 changes: 2 additions & 4 deletions R/globalVars.R
Original file line number Diff line number Diff line change
@@ -1,9 +1,7 @@
utils::globalVariables(c(
"Period", "Period1", "Period2", "Run",
"Scenario", "Year", "bands", "conn",
"fullnm", "laynum", "list_gcm_ts",
"list_period", "list_years_ts", "mod", "nameque",
"numlay", "period", "projID", "rastque", "run", "sprc", "var"
"Scenario", "Year", "fullnm", "laynum",
"mod", "numlay", "period", "run", "var"
))

#
Loading
Loading