[R] Problem with combining monthly nc files into a yearly file (era5 climate data)
Roy Mendelssohn - NOAA Federal
roy@mende|@@ohn @end|ng |rom no@@@gov
Fri Jun 21 16:14:37 CEST 2024
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 using 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 using 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.
More information about the R-help
mailing list