Dear R-help List,
I am currently trying to run a code which is available on Zenodo
(https://zenodo.org/records/10997880 - 02_MicroClimModel.R).
The code downloads yearly era5 climate data. Unfortunately, the limit to
download these nc-files was recently reduced to 60000. Therefore, I can not
download the yearly file anymore. I have solved this by rewriting the code, so
that it downloads 12 monthly files.
However, I have not been able to combine these 12 monthly nc-files into one
yearly file. The code gives me errors if I continue running it. I assume that
the combination was not successful and might have messed up the format. I would
greatly appreciate any advice on how to convert these monthly nc-files into one
yearly file.
Thank you very much in advance!
Here is the full code:
' *****************************************************************
#' ~~ STEP 01 DOWNLOADING & PROCESSING HOURLY CLIMATE DATA
# Install the remotes package if not already installed
if (!requireNamespace("remotes", quietly = TRUE)) {
install.packages("remotes")
}
# Install packages from CRAN
install.packages(c("terra", "raster", "ncdf4", "lubridate"))
install.packages("lutz")
#install dependencies for microclima
remotes::install_github("ropensci/rnoaa")
# Install packages from GitHub
remotes::install_github("dklinges9/mcera5")
remotes::install_github("ilyamaclean/microclima")
remotes::install_github("ilyamaclean/microclimf")
#' ~~ Required libraries:
require(terra)
require(raster)
require(mcera5) # https://github.com/dklinges9/mcera5
require(ncdf4)
require(microclima) # https://github.com/ilyamaclean/microclima
require(microclimf) # https://github.com/ilyamaclean/microclimf
require(ecmwfr)
require(lutz)
require(lubridate)
# Set paths and year of interest
pathtodata <- "F:/Dat/"
pathtoera5 <- paste0(pathtodata, "era5/")
year <- 2019
# Set user credentials for CDS API (you have to first register and insert here
your UID and API key at https://cds.climate.copernicus.eu/user/register and
allow downloads)
uid <- "xxx"
cds_api_key <- "xxx"
ecmwfr::wf_set_key(user = uid, key = cds_api_key, service = "cds")
# Define the spatial extent for your tile
xmn <- 18.125
xmx <- 22.875
ymn <- -1.625
ymx <- 1.875
#HERE STARTS THE SECTION WHERE I AM DOWNLOADING MONTHLY FILES
# Define the temporal extent of the run
start_time <- lubridate::ymd(paste0(year, "-01-01"))
end_time <- lubridate::ymd(paste0(year, "-12-31"))
# Function to build and send requests for each month
request_era5_monthly <- function(year, month, uid, xmn, xmx, ymn, ymx,
out_path) {
# Define the start and end times for the month
st_time <- lubridate::ymd(paste0(year, "-", sprintf("%02d", month), "-01"))
en_time <- st_time + months(1) - days(1)
# Create the file prefix and request
file_prefix <- paste0("era5_reanalysis_", year, "_", sprintf("%02d", month))
req <- build_era5_request(xmin = xmn, xmax = xmx, ymin = ymn, ymax = ymx,
start_time = st_time, end_time = en_time,
outfile_name = file_prefix)
# Send the request and save the data
request_era5(request = req, uid = uid, out_path = out_path, overwrite = TRUE)
}
# Loop over each month and request data
for (month in 1:12) {
request_era5_monthly(year, month, uid, xmn, xmx, ymn, ymx, pathtoera5)
}
#HERE I AM EXPLORING ONE EXEMPLARY MONTHLY NC FILE
file_path <- paste0(pathtoera5, "era5_reanalysis_2019_01_2019.nc")
nc <- nc_open(file_path)
# List all variables
print(nc)
# List all variable names in the NetCDF file
var_names <- names(nc$var)
print(var_names)
checkJan <- raster(paste0(pathtoera5, "era5_reanalysis_2019_01_2019.nc"))
print(checkJan)
opencheckJan <- getValues(checkJan)
opencheckJan
#HERE IS THE PROBLEM, I AM TRYING TO COMBINE THESE MONTHL NC FILES
combine_era5_yearly <- function(year, pathtoera5, outfile) {
# List of monthly files
monthly_files <- list.files(pathtoera5, pattern = paste0("era5_reanalysis_",
year, "_\\d{2}_", year, "\\.nc"), full.names = TRUE)
if (length(monthly_files) == 0) {
stop("No monthly files found")
}
# Initialize lists to store data
lons <- NULL
lats <- NULL
time <- NULL
t2m <- list()
d2m <- list()
sp <- list()
u10 <- list()
v10 <- list()
tp <- list()
tcc <- list()
msnlwrf <- list()
msdwlwrf <- list()
fdir <- list()
ssrd <- list()
lsm <- list()
# Read each monthly file and extract variables
for (file in monthly_files) {
nc <- nc_open(file)
if (is.null(lons)) {
lons <- ncvar_get(nc, "longitude")
lats <- ncvar_get(nc, "latitude")
time <- ncvar_get(nc, "time")
} else {
time <- c(time, ncvar_get(nc, "time"))
}
t2m <- c(t2m, list(ncvar_get(nc, "t2m")))
d2m <- c(d2m, list(ncvar_get(nc, "d2m")))
sp <- c(sp, list(ncvar_get(nc, "sp")))
u10 <- c(u10, list(ncvar_get(nc, "u10")))
v10 <- c(v10, list(ncvar_get(nc, "v10")))
tp <- c(tp, list(ncvar_get(nc, "tp")))
tcc <- c(tcc, list(ncvar_get(nc, "tcc")))
msnlwrf <- c(msnlwrf, list(ncvar_get(nc, "msnlwrf")))
msdwlwrf <- c(msdwlwrf, list(ncvar_get(nc, "msdwlwrf")))
fdir <- c(fdir, list(ncvar_get(nc, "fdir")))
ssrd <- c(ssrd, list(ncvar_get(nc, "ssrd")))
lsm <- c(lsm, list(ncvar_get(nc, "lsm")))
nc_close(nc)
}
# Combine the data for each variable
t2m <- do.call(c, t2m)
d2m <- do.call(c, d2m)
sp <- do.call(c, sp)
u10 <- do.call(c, u10)
v10 <- do.call(c, v10)
tp <- do.call(c, tp)
tcc <- do.call(c, tcc)
msnlwrf <- do.call(c, msnlwrf)
msdwlwrf <- do.call(c, msdwlwrf)
fdir <- do.call(c, fdir)
ssrd <- do.call(c, ssrd)
lsm <- do.call(c, lsm)
# Create a new NetCDF file for the entire year
outfile <- paste0(pathtoera5, "era5_reanalysis_", year, ".nc")
dim_lon <- ncdim_def("longitude", "degrees_east", lons)
dim_lat <- ncdim_def("latitude", "degrees_north", lats)
dim_time <- ncdim_def("time", "hours since 1900-01-01 00:00:00", time,
unlim=TRUE)
# Define variables
var_t2m <- ncvar_def("t2m", "K", list(dim_lon, dim_lat, dim_time), -9999)
var_d2m <- ncvar_def("d2m", "K", list(dim_lon, dim_lat, dim_time), -9999)
var_sp <- ncvar_def("sp", "Pa", list(dim_lon, dim_lat, dim_time), -9999)
var_u10 <- ncvar_def("u10", "m/s", list(dim_lon, dim_lat, dim_time), -9999)
var_v10 <- ncvar_def("v10", "m/s", list(dim_lon, dim_lat, dim_time), -9999)
var_tp <- ncvar_def("tp", "m", list(dim_lon, dim_lat, dim_time), -9999)
var_tcc <- ncvar_def("tcc", "1", list(dim_lon, dim_lat, dim_time), -9999)
var_msnlwrf <- ncvar_def("msnlwrf", "W/m^2", list(dim_lon, dim_lat,
dim_time), -9999)
var_msdwlwrf <- ncvar_def("msdwlwrf", "W/m^2", list(dim_lon, dim_lat,
dim_time), -9999)
var_fdir <- ncvar_def("fdir", "J/m^2", list(dim_lon, dim_lat, dim_time),
-9999)
var_ssrd <- ncvar_def("ssrd", "J/m^2", list(dim_lon, dim_lat, dim_time),
-9999)
var_lsm <- ncvar_def("lsm", "1", list(dim_lon, dim_lat, dim_time), -9999)
# Create the file
ncout <- nc_create(outfile, list(var_t2m, var_d2m, var_sp, var_u10, var_v10,
var_tp, var_tcc, var_msnlwrf, var_msdwlwrf, var_fdir, var_ssrd, var_lsm))
# Write data to the new file
ncvar_put(ncout, var_t2m, t2m)
ncvar_put(ncout, var_d2m, d2m)
ncvar_put(ncout, var_sp, sp)
ncvar_put(ncout, var_u10, u10)
ncvar_put(ncout, var_v10, v10)
ncvar_put(ncout, var_tp, tp)
ncvar_put(ncout, var_tcc, tcc)
ncvar_put(ncout, var_msnlwrf, msnlwrf)
ncvar_put(ncout, var_msdwlwrf, msdwlwrf)
ncvar_put(ncout, var_fdir, fdir)
ncvar_put(ncout, var_ssrd, ssrd)
ncvar_put(ncout, var_lsm, lsm)
# Define and write longitude and latitude variables
ncvar_put(ncout, "longitude", lons)
ncvar_put(ncout, "latitude", lats)
ncatt_put(ncout, "longitude", "units", "degrees_east")
ncatt_put(ncout, "latitude", "units", "degrees_north")
# Define and write time variable
ncvar_put(ncout, "time", time)
ncatt_put(ncout, "time", "units", "hours since 1900-01-01 00:00:00")
ncatt_put(ncout, "time", "calendar", "gregorian")
# Global attributes
ncatt_put(ncout, 0, "title", paste0("ERA5 reanalysis data for ", year))
ncatt_put(ncout, 0, "source", "ECMWF ERA5")
# Close the NetCDF file
nc_close(ncout)
}
# Example usage:
outfile <- paste0(pathtoera5, "era5_reanalysis_", year, ".nc")
combine_era5_yearly(year, pathtoera5, outfile)
#HERE IS THE REST OF THE CODE WHICH REQUIRES THE YEARLY FILE
#' Process Hourly Climate Data >>>
file <- paste0(pathtoera5,"era5_reanalysis_",year,".nc")
clim <- nc_open(file)
#' create a template to crop input dataset to for step 2: '02_VegParms.R'
test <- raster::brick(file, varname = "t2m")
t_array <- as.array(test[[1]])
ext_r <- ext(raster::extent(test))
r <- rast(t_array, crs = "EPSG:4326", ext = ext_r)
#' Get coordinates & time:
lons <- ncdf4::ncvar_get(clim, varid = "longitude")
lats <- ncdf4::ncvar_get(clim, varid = "latitude")
time <- ncdf4::ncvar_get(clim, "time")
x_dim <- length(lons)
y_dim <- length(lats)
z_dim <- length(time)
#' Assign a local timezone:
tmz <- lutz::tz_lookup_coords(lats[length(lats)/2], lons[length(lons)/2],
method = 'fast')
origin <- as.POSIXlt("1900-01-01 00:00:00", tz = "UTC")
UTC_tme <- origin + as.difftime(time, units = "hours")
UTC_tme <- as.POSIXlt(UTC_tme, tz = "UTC")
local_tme <- lubridate::with_tz(UTC_tme, tzone = tmz)
jd <- microctools::jday(tme = UTC_tme)
lt <- local_tme$hour + local_tme$min/60 + local_tme$sec/3600
#' Create empty climate variable arrays:
#' These are 3-D arrays with time (hours) in the 3rd dimension.
t_a <- array(data = NA, c(y_dim, x_dim, z_dim))
t_sh <- array(data = NA, c(y_dim, x_dim, z_dim))
t_pa <- array(data = NA, c(y_dim, x_dim, z_dim))
t_ws <- array(data = NA, c(y_dim, x_dim, z_dim))
t_wd <- array(data = NA, c(y_dim, x_dim, z_dim))
t_se <- array(data = NA, c(y_dim, x_dim, z_dim))
t_nl <- array(data = NA, c(y_dim, x_dim, z_dim))
t_ul <- array(data = NA, c(y_dim, x_dim, z_dim))
t_dl <- array(data = NA, c(y_dim, x_dim, z_dim))
t_rd <- array(data = NA, c(y_dim, x_dim, z_dim))
t_rdf <- array(data = NA, c(y_dim, x_dim, z_dim))
t_sz <- array(data = NA, c(y_dim, x_dim, z_dim))
p_a <- array(data = NA, c(y_dim, x_dim, length(local_tme)/24)) # note: rainfall
recorded daily.
#' Fill empty arrays with processed era5 data:
#' Use the 'extract_clim' function from the mcera5 package to convert era5 data
to microclimf-ready data.
for(i in 1:y_dim){ # for each row in the new array.
for(j in 1:x_dim){ # for each column in the new array.
long <- lons[j]
lat <- lats[i]
climate <- extract_clim(file, long, lat, start_time = UTC_tme[1], end_time
= UTC_tme[length(UTC_tme)])
t_a[i,j,] <- climate$temperature
t_sh[i,j,] <- climate$humidity
t_pa[i,j,] <- climate$pressure
t_ws[i,j,] <- climate$windspeed
t_wd[i,j,] <- climate$winddir
t_se[i,j,] <- climate$emissivity
t_nl[i,j,] <- climate$netlong
t_ul[i,j,] <- climate$uplong
t_dl[i,j,] <- climate$downlong
t_rd[i,j,] <- climate$rad_dni
t_rdf[i,j,] <- climate$rad_dif
t_sz[i,j,] <- climate$szenith
} # end column j
} # end row i
#' Repeat for daily rainfall:
#' Use the 'extract_precip' function from the mcera5 package to convert era5
data to microclimf-ready data.
for(i in 1:y_dim){ # for each row in the new array.
for(j in 1:x_dim){ # for each column in the new array.
long <- lons[j]
lat <- lats[i]
precip <- extract_precip(file, long, lat, start_time = UTC_tme[1], end_time
= UTC_tme[length(UTC_tme)])
p_a[i,j,] <- precip
} # end column j
} # end row i
#' Additional processing for microclimf inputs:
#' Calculate Solar Index...
si <- array(data = NA, dim = c(y_dim, x_dim, z_dim))
for(a in 1:nrow(si)){
for(b in 1:ncol(si)){
x <- lons[b]
y <- lats[a]
s = microclima::siflat(lubridate::hour(local_tme), y, x, jd)
si[a,b,] <- s
}
}
#' Calculate Global Horizontal Irradiance (GHI) from Direct Normal Irradiance
(DNI)
#' and convert units from MJh/m^2 to kWh/m^2
raddr <- (t_rd * si)/0.0036
difrad <- t_rdf/0.0036
#' Cap diffuse radiation data (Cannot be less than 0)
difrad[difrad < 0] <- 0
#' Calculate shortwave radiation:
#' Sum Global Horizontal Irradiance (GHI) and Diffuse Radiation.
swrad <- raddr + difrad
#' Cap shortwave radiation between 0 > sw < 1350 (lower than the solar constant)
swrad[swrad < 0] <- 0
swrad[swrad > 1350] <- 1350
#' Calculate relative humidity
#' Using specific humidity, temperature and pressure.
t_rh <- array(data = NA, c(y_dim, x_dim, z_dim))
for(i in 1:nrow(t_rh)){
for(j in 1:ncol(t_rh)){
rh <- microclima::humidityconvert(t_sh[i,j,],intype = "specific", tc =
t_a[i,j,], p = t_pa[i,j,])
rh <- rh$relative
rh[rh > 100] <- 100
t_rh[i,j,] <- rh
}
}
#' Convert pressure untis from Pa to kPa:
t_pr <- t_pa/1000
#' Create final climate data set to drive microclimate model:
#' Note: keep the nomenclature as shown here for microclimf, see
microclimf::climdat for example names.
climdat <- list(tme = local_tme, obs_time = UTC_tme,
temp = t_a, relhum = t_rh,
pres = t_pr, swrad = swrad,
difrad = difrad, skyem = t_se,
windspeed = t_ws, winddir = t_wd)
#' Save data:
pathout <- "F:/Dat/era5/"
saveRDS(climdat, paste0(pathout,"climdat_",year,".RDS"))
saveRDS(p_a, paste0(pathout,"rainfall_",year,".RDS"))
tile_no <- "01"
writeRaster(r, paste0(pathout,"tile_",tile_no,".tif"))
#HERE IS ADDITIONAL INFORMATION ON ONE MONTHLY NC FILE:
12 variables (excluding dimension variables):
short t2m[longitude,latitude,time]
scale_factor: 0.000250859493618673
add_offset: 301.508114316347
_FillValue: -32767
missing_value: -32767
units: K
long_name: 2 metre temperature
short d2m[longitude,latitude,time]
scale_factor: 0.000189842033307647
add_offset: 296.056545703983
_FillValue: -32767
missing_value: -32767
units: K
long_name: 2 metre dewpoint temperature
short sp[longitude,latitude,time]
scale_factor: 0.0470135275357454
add_offset: 96477.3202432362
_FillValue: -32767
missing_value: -32767
units: Pa
long_name: Surface pressure
standard_name: surface_air_pressure
short u10[longitude,latitude,time]
scale_factor: 0.000152449582891444
add_offset: 0.590744087708554
_FillValue: -32767
missing_value: -32767
units: m s**-1
long_name: 10 metre U wind component
short v10[longitude,latitude,time]
scale_factor: 0.00013693746249206
add_offset: 0.66616840871016
_FillValue: -32767
missing_value: -32767
units: m s**-1
long_name: 10 metre V wind component
short tp[longitude,latitude,time]
scale_factor: 2.85070516901134e-07
add_offset: 0.00934062055678257
_FillValue: -32767
missing_value: -32767
units: m
long_name: Total precipitation
short tcc[longitude,latitude,time]
scale_factor: 1.52594875864068e-05
add_offset: 0.499992370256207
_FillValue: -32767
missing_value: -32767
units: (0 - 1)
long_name: Total cloud cover
standard_name: cloud_area_fraction
short msnlwrf[longitude,latitude,time]
scale_factor: 0.00173717121168915
add_offset: -56.7456195621683
_FillValue: -32767
missing_value: -32767
units: W m**-2
long_name: Mean surface net long-wave radiation flux
short msdwlwrf[longitude,latitude,time]
scale_factor: 0.0012878582820392
add_offset: 410.789761344296
_FillValue: -32767
missing_value: -32767
units: W m**-2
long_name: Mean surface downward long-wave radiation flux
short fdir[longitude,latitude,time]
scale_factor: 46.9767598004059
add_offset: 1539240.5116201
_FillValue: -32767
missing_value: -32767
units: J m**-2
long_name: Total sky direct solar radiation at surface
short ssrd[longitude,latitude,time]
scale_factor: 54.2183022294111
add_offset: 1776516.89084889
_FillValue: -32767
missing_value: -32767
units: J m**-2
long_name: Surface short-wave (solar) radiation downwards
standard_name: surface_downwelling_shortwave_flux_in_air
short lsm[longitude,latitude,time]
scale_factor: 9.55416624213488e-06
add_offset: 0.686938634743966
_FillValue: -32767
missing_value: -32767
units: (0 - 1)
long_name: Land-sea mask
standard_name: land_binary_mask
3 dimensions:
longitude Size:21
units: degrees_east
long_name: longitude
latitude Size:16
units: degrees_north
long_name: latitude
time Size:744
units: hours since 1900-01-01 00:00:00.0
long_name: time
calendar: gregorian
2 global attributes:...
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
