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@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@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. ______________________________________________ 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.
