From f5f8711709afb2489b7c71f60786ebcf90894f3d Mon Sep 17 00:00:00 2001 From: Ghislain Vieilledent Date: Sun, 10 Sep 2023 22:02:39 +0200 Subject: [PATCH] Replacing stars by terra in get_chelsa_current.R --- R/get_chelsa_current.R | 146 +++++++++-------------- R/get_nodata_dtype_chelsa.R | 25 ++-- data/nodata_chelsa.rda | Bin 377 -> 377 bytes tests/testthat/test-get_chelsa_current.R | 3 +- 4 files changed, 73 insertions(+), 101 deletions(-) diff --git a/R/get_chelsa_current.R b/R/get_chelsa_current.R index b312134..de7acc4 100644 --- a/R/get_chelsa_current.R +++ b/R/get_chelsa_current.R @@ -111,106 +111,63 @@ get_chelsa_current <- function(extent_latlon, extent_proj, EPSG_proj, destinatio "envicloud/chelsa/chelsa_V2/GLOBAL/", "climatologies/1981-2010") - progress_bar <- 0 - nb_var_download <- 12 * 6 - cat("Downloading tasmin, tasmax, tas, pr, clt, and pet_penman\n") - pb = txtProgressBar(min = 0, max = nb_var_download, initial = 0) - for (m in stringr::str_pad(1:12, width = 2, pad = "0")) { - - ## Monthly minimum temperature (°C). - ifile <- glue::glue("/vsicurl/{url_base_chelsa}/tasmin/", - "CHELSA_tasmin_{m}_1981-2010_V.2.1.tif") - ofile <- file.path(destination, 'data_raw', 'chelsa_v2_1', 'temp', - paste0('tasmin', m, '.tif')) - gdal_utils_translate(ifile, ofile, extent_gdal_translate, opts=c("-ot", "Int32")) - progress_bar <- progress_bar + 1 - setTxtProgressBar(pb, progress_bar) - - ## Monthly maximum temperature (°C). - ifile <- glue::glue("/vsicurl/{url_base_chelsa}/tasmax/", - "CHELSA_tasmax_{m}_1981-2010_V.2.1.tif") - ofile <- file.path(destination, 'data_raw', 'chelsa_v2_1', 'temp', - paste0('tasmax', m, '.tif')) - gdal_utils_translate(ifile, ofile, extent_gdal_translate, opts=c("-ot", "Int32")) - progress_bar <- progress_bar + 1 - setTxtProgressBar(pb, progress_bar) - - ## Monthly average temperature (°C). - ifile <- glue::glue("/vsicurl/{url_base_chelsa}/tas/", - "CHELSA_tas_{m}_1981-2010_V.2.1.tif") - ofile <- file.path(destination, 'data_raw', 'chelsa_v2_1', 'temp', - paste0('tas', m, '.tif')) - gdal_utils_translate(ifile, ofile, extent_gdal_translate, opts=c("-ot", "Int32")) - progress_bar <- progress_bar + 1 - setTxtProgressBar(pb, progress_bar) - - ## Monthly precipitation (mm ~ kg/m2). - ifile <- glue::glue("/vsicurl/{url_base_chelsa}/pr/", - "CHELSA_pr_{m}_1981-2010_V.2.1.tif") - ofile <- file.path(destination, 'data_raw', 'chelsa_v2_1', 'temp', - paste0('pr', m, '.tif')) - gdal_utils_translate(ifile, ofile, extent_gdal_translate, opts=c("-ot", "Int32")) - progress_bar <- progress_bar + 1 - setTxtProgressBar(pb, progress_bar) - - ## Monthly cloud area fraction - ifile <- glue::glue("/vsicurl/{url_base_chelsa}/clt/", - "CHELSA_clt_{m}_1981-2010_V.2.1.tif") - ofile <- file.path(destination, 'data_raw', 'chelsa_v2_1', 'temp', - paste0('clt', m, '.tif')) - gdal_utils_translate(ifile, ofile, extent_gdal_translate, opts=c("-ot", "Int32")) - progress_bar <- progress_bar + 1 - setTxtProgressBar(pb, progress_bar) - - ## Monthly pet_penman - # https://www.fao.org/3/x0490e/x0490e06.htm - ifile <- glue::glue("/vsicurl/{url_base_chelsa}/pet/", - "CHELSA_pet_penman_{m}_1981-2010_V.2.1.tif") - ofile <- file.path(destination, 'data_raw', 'chelsa_v2_1', 'temp', - paste0('pet_penman', m, '.tif')) - gdal_utils_translate(ifile, ofile, extent_gdal_translate, opts=c("-ot", "Int32")) - progress_bar <- progress_bar + 1 - setTxtProgressBar(pb, progress_bar) + ## Climatic variables + var <- c("clt", "pet", "pr", "tas", "tasmax", "tasmin") + varr <- c("clt", "pet_penman", "pr", "tas", "tasmax", "tasmin") + nvar <- length(var) + + iter_bar <- 0 + nb_var_download <- nvar * 12 + cat("Downloading clt, pet_penman, pr, tas, tasmax, tasmin", "\n") + pb = txtProgressBar(min = 1, max = nb_var_download, width=30, style=3) + ## Loop on variables + for (v in 1:nvar) { + ## Loop on months + for (m in stringr::str_pad(1:12, width = 2, pad = "0")) { + ifile <- glue::glue("/vsicurl/{url_base_chelsa}/{var[v]}/", + "CHELSA_{varr[v]}_{m}_1981-2010_V.2.1.tif") + ofile <- file.path(destination, "data_raw", "chelsa_v2_1", "temp", + paste0(varr[v], m, ".tif")) + gdal_utils_translate(ifile, ofile, extent_gdal_translate, + opts=c("-ot", "Int32")) + iter_bar <- iter_bar + 1 + setTxtProgressBar(pb, iter_bar) + } } close(pb) ## Bioclimatic variables - ## See https://chelsa-climate.org/wp-admin/download-page/CHELSA_tech_specification_V2.pdf for details - cat("Downloading bioclimatic variables\n") + ## https://chelsa-climate.org/wp-admin/download-page/CHELSA_tech_specification_V2.pdf + cat("Downloading bioclimatic variables", "\n") + pb = txtProgressBar(min = 1, max = 19, width=30, style=3) for (i in 1:19) { ifile <- glue::glue("/vsicurl/{url_base_chelsa}/bio/", "CHELSA_bio{i}_1981-2010_V.2.1.tif") ofile <- file.path(destination, 'data_raw', 'chelsa_v2_1', 'temp', paste0('bio', stringr::str_pad(i, width = 2, pad = '0'), '.tif')) # Data type and nodata - # There are errors in Chelsa metadata as NoData and type do not correspond - # eg: UInt16 and -99999 - metadata <- sf::gdal_utils("gdalinfo", ifile, quiet=TRUE) - # data_type <- regmatches(metadata, regexpr("Type=[[:graph:]]+", metadata)) - # data_type <- sub("Type=", "", data_type) - nodata_val <- regmatches(metadata, regexpr("NoData[[:space:]]Value=[[:graph:]]+", metadata)) - nodata_val <- as.numeric(sub("NoData Value=", "", nodata_val)) - # cat(data_type, nodata_val, "\n") - if (nodata_val == 65535) {dtype <- "UInt16"} - if (nodata_val == -99999) {dtype <- "Int32"} - if (nodata_val > 3.4e+38) {dtype <- "Float32"} + if (i %in% c(8:11, 16:19)) {dtype <- "UInt16"} # nodata_val == 65535 + if (i %in% c(1, 2, 4:7, 12:14)) {dtype <- "Int32"} # nodata_val == -99999 + if (i %in% c(3, 15)) {dtype <- "Float32"} # nodata_val == 3.4028234663852886e+38 # Download gdal_utils_translate(ifile, ofile, extent_gdal_translate, opts=c("-ot", dtype)) + setTxtProgressBar(pb, i) } + close(pb) ## ================================ - ## Reproject and stack per variable + ## Reproject ## ================================ - cat("Reprojecting rasters and stacking monthly rasters per variable\n") - for (var in c("tasmin", "tasmax", "tas", "pr", "bio", "clt", "pet_penman")) { - + cat("Reprojecting rasters", "\n") + + var_list <- c("tas", "tasmin", "tasmax", "pr", "pet_penman", "clt", "bio") + for (var in var_list) { idir <- file.path(destination, "data_raw", "chelsa_v2_1", "temp") - tif_files <- list.files(idir, pattern = paste0(var, "[0-9]{2}\\.tif"), full.names = TRUE) - + tif_files <- list.files(idir, pattern = paste0(var, "[0-9]{2}\\.tif"), + full.names = TRUE) for (i in 1:length(tif_files)) { - ifile <- tif_files[i] ofile <- gsub(".tif", "_res.tif", tif_files[i]) opts <- glue("-tr {resol} {resol} -te {extent_proj_string} ", @@ -221,16 +178,24 @@ get_chelsa_current <- function(extent_latlon, extent_proj, EPSG_proj, destinatio options=unlist(strsplit(opts, " ")), quiet=TRUE) } + } + + ## ================================ + ## Stack per variable + ## ================================ - # Stack per variable + cat("Stack per variable", "\n") + for (var in var_list) { idir <- file.path(destination, "data_raw", "chelsa_v2_1", "temp") - tif_files_res <- list.files(idir, pattern = paste0(var, "[0-9]{2}_res\\.tif"), full.names = TRUE) + tif_files_res <- list.files(idir, pattern = paste0(var, "[0-9]{2}_res\\.tif"), + full.names = TRUE) r <- terra::rast(sort(tif_files_res)) names(r) <- paste0(var, 1:length(tif_files_res)) # Export ofile <- file.path(destination, "data_raw", "chelsa_v2_1", paste0(var, "_res.tif")) - terra::writeRaster(r, gdal = c("COMPRESS=LZW","PREDICTOR=2"), progress = 0, overwrite = TRUE, + terra::writeRaster(r, gdal = c("COMPRESS=LZW","PREDICTOR=2"), + progress = 0, overwrite = TRUE, datatype = "INT2S", filename = ofile) } @@ -238,16 +203,15 @@ get_chelsa_current <- function(extent_latlon, extent_proj, EPSG_proj, destinatio ## Create raster stack ## ============================== - cat("Create raster stack of monthly variables and bioclimatic variables\n") + cat("Create raster stack of monthly variables and bioclimatic variables", "\n") # Stack tasmin, tasmax, tas, pr, clt, pet_penman, and bio - files.tif <- file.path(destination, "data_raw", "chelsa_v2_1", - paste0(c("tasmin", "tasmax", "tas", "pr", "clt", "pet_penman", "bio"), "_res.tif")) - r <- c(read_stars(files.tif[1]), read_stars(files.tif[2]), read_stars(files.tif[3]), read_stars(files.tif[4]), - read_stars(files.tif[5]), read_stars(files.tif[6]), read_stars(files.tif[7]), along = "band") + tif_files <- file.path(destination, "data_raw", "chelsa_v2_1", + paste0(var_list, "_res.tif")) + r <- terra::rast(tif_files) ofile <- file.path(destination, "data_raw", "chelsa_v2_1", "clim_res.tif") - write_stars(obj = r, options = c("COMPRESS=LZW", "PREDICTOR=2"), NA_value = nodat, - type = "Int16", dsn = ofile) - rm(r) + terra::writeRaster(r, gdal = c("COMPRESS=LZW","PREDICTOR=2"), + progress = 0, overwrite = TRUE, + datatype = "INT2S", filename = ofile) ## =================================================== ## Compute water deficit (cwd and ndw) with Penman ETP diff --git a/R/get_nodata_dtype_chelsa.R b/R/get_nodata_dtype_chelsa.R index 51a2b9d..244b265 100644 --- a/R/get_nodata_dtype_chelsa.R +++ b/R/get_nodata_dtype_chelsa.R @@ -46,6 +46,7 @@ get_nodata_dtype_chelsa <- function (nmonths=1, biovar=TRUE, verbose=TRUE) { # Climatic variables var <- c("clt", "pet", "pr", "tas", "tasmax", "tasmin", bio) varr <- c("clt", "pet_penman", "pr", "tas", "tasmax", "tasmin", bioo) + nvar <- length(var) # Base url for Chelsa url_base_chelsa <- paste0( @@ -55,9 +56,9 @@ get_nodata_dtype_chelsa <- function (nmonths=1, biovar=TRUE, verbose=TRUE) { # Dataframe to store results df <- data.frame() - + # Loop on variables - for (v in 1:length(var)) { + for (v in 1:nvar) { # Message if (verbose) { cat(glue("Get nodata value and data type for: \"{varr[v]}\"."), "\n") @@ -69,13 +70,15 @@ get_nodata_dtype_chelsa <- function (nmonths=1, biovar=TRUE, verbose=TRUE) { ifile <- glue::glue("/vsicurl/{url_base_chelsa}/{var[v]}/", "CHELSA_{varr[v]}_{m}_1981-2010_V.2.1.tif") # Get metadata using gdalinfo and regular expressions - metadata <- sf::gdal_utils("gdalinfo", ifile, quiet=TRUE) - data_type <- regmatches(metadata, regexpr("Type=[[:graph:]]+", metadata)) + metadata <- sf::gdal_utils("gdalinfo", ifile, + config_options=c(GTIFF_SRS_SOURCE="EPSG"), + quiet=TRUE) + data_type <- regmatches(metadata, regexpr("Type=[[:alnum:]]+", metadata)) data_type <- sub("Type=", "", data_type) reg_expr <- regexpr("NoData[[:space:]]Value=[[:graph:]]+", metadata) nodata_val <- regmatches(metadata, reg_expr) - nodata_val <- as.numeric(sub("NoData Value=", "", nodata_val)) - df <- rbind(df, c(varr[v], as.numeric(m), data_type, nodata_val)) + nodata_val <- sub("NoData Value=", "", nodata_val) + df <- rbind(df, c(varr[v], m, data_type, nodata_val)) } } if (varr[v] %in% bioo) { @@ -83,17 +86,21 @@ get_nodata_dtype_chelsa <- function (nmonths=1, biovar=TRUE, verbose=TRUE) { ifile <- glue::glue("/vsicurl/{url_base_chelsa}/{var[v]}/", "CHELSA_{varr[v]}_1981-2010_V.2.1.tif") # Get metadata using gdalinfo and regular expressions - metadata <- sf::gdal_utils("gdalinfo", ifile, quiet=TRUE) - data_type <- regmatches(metadata, regexpr("Type=[[:graph:]]+", metadata)) + metadata <- sf::gdal_utils("gdalinfo", ifile, + config_options=c(GTIFF_SRS_SOURCE="EPSG"), + quiet=TRUE) + data_type <- regmatches(metadata, regexpr("Type=[[:alnum:]]+", metadata)) data_type <- sub("Type=", "", data_type) reg_expr <- regexpr("NoData[[:space:]]Value=[[:graph:]]+", metadata) nodata_val <- regmatches(metadata, reg_expr) - nodata_val <- as.numeric(sub("NoData Value=", "", nodata_val)) + nodata_val <- sub("NoData Value=", "", nodata_val) df <- rbind(df, c(varr[v], NA, data_type, nodata_val)) } } # Rename columns names(df) <- c("var", "month", "dtype", "nodata") + df$m <- as.numeric(df$m) + df$nodata <- as.numeric(df$nodata) # Return dataframe with results return(df) diff --git a/data/nodata_chelsa.rda b/data/nodata_chelsa.rda index 07d3054e84b50224db3f064b7713a8dc5db06166..575c0ec7341e3ac02af084f5f195eeadf482a032 100644 GIT binary patch literal 377 zcmV-<0fzn`iwFP!000001I?7pO2a@9$EQit(m`uW=m{gZ%y9L(=ivFC-1TaABO{_hB2M;!CvGJO(%u+;ei z?_joBXDc@m^Q0ax)AafDaizvpRo78F*q9;}XC&Txd5x@>hU1nATD-_p7L{vI=H^0mpe zM0Bx^*g_DJ5QI0+5IRWQv)6`D3K$5>WBF)_Bq0TM)-lORY9 z*&+b!V4S1@>;<8%24hYY)?*d$(QFxn1i~*;z{X`Nt4ED3^ayo@g?0LuaD)d)MA}fw z2$+dlw1G%82;`jS07`J~4+lQPg+dRFpCBu?fVS?kt66#6UxU{gI4JJWA%Ojg=_hAN zg_^fH<96b`_oNYU*%nH(CdC#Ibh|w`fT@)r(E@!uv}fN{9I&G!K4wiY3MVqMk^OG! Xj=8V=86hT_ApaL~ML1B9UDRjQWz(j+ diff --git a/tests/testthat/test-get_chelsa_current.R b/tests/testthat/test-get_chelsa_current.R index 0c6c32f..708afb8 100644 --- a/tests/testthat/test-get_chelsa_current.R +++ b/tests/testthat/test-get_chelsa_current.R @@ -2,7 +2,8 @@ library(gecevar) iso <- "REU" epsg <- 32740 r <- get_aoi_extent(EPSG_proj=epsg, - country_iso=iso) + country_iso=iso, + rm_output_dir=FALSE) extent_latlon <- r$extent_latlon extent_proj <- r$extent_proj resol <- 1000