R help - Jun 2024 - Problem with combining monthly nc files into a yearly file (era5 climate data)

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:...

Hi Leni:

You forget to post the important part - the errors you have been getting and if
you have the errors isolated to particular lines in the code.

HTH,

-Roy

> On Jun 21, 2024, at 3:59?AM, Leni Koehnen via R-help <r-help at
r-project.org> wrote:
> 
> 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 at 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.