From b5b274a4f84dcf677d752a31d0d11a34577a9dbe Mon Sep 17 00:00:00 2001 From: nperez Date: Thu, 27 Mar 2025 10:57:40 +0100 Subject: [PATCH 1/2] Improvements to orography correction --- modules/Loading/R/orography_correction.R | 123 ++++++++++++++--------- 1 file changed, 77 insertions(+), 46 deletions(-) diff --git a/modules/Loading/R/orography_correction.R b/modules/Loading/R/orography_correction.R index 3e12c4c2..aa1ae035 100644 --- a/modules/Loading/R/orography_correction.R +++ b/modules/Loading/R/orography_correction.R @@ -1,75 +1,106 @@ -## This function is to apply an orographic temperature correction between the -## exp and obs datasets. This correction is relevant for the mean_bias metric. - +## TODO: remove paths to personal scratchs source("modules/Loading/R/get_regrid_params.R") source("modules/Loading/R/check_latlon.R") -orography_correction <- function(recipe, data) { - +# This correction can only be applied to temperature. +# This correction is relevant for the mean_bias metric. +# correction_constant is negative +# -6.5ºK/km(ºC/km for the environmental lapse rate correction +# -9.8ºK/km (�C/km) for the adiabatic correction +orography_correction <- function(recipe, data, correction_constant = -6.5) { + ref.name <- recipe$Analysis$Datasets$Reference$name exp.name <- recipe$Analysis$Datasets$System$name lats.min <- recipe$Analysis$Region$latmin lats.max <- recipe$Analysis$Region$latmax lons.min <- recipe$Analysis$Region$lonmin lons.max <- recipe$Analysis$Region$lonmax - - archive <- get_archive(recipe) + + archive <- read_yaml("conf/archive_seasonal.yml")[[recipe$Run$filesystem]] # Define regrid parameters: regrid_params <- get_regrid_params(recipe, archive) - + circularsort <- check_latlon(lats.min, lats.max, lons.min, lons.max) - + ## Load exp orography - orography_exp <- Start(dat = archive$System[[exp.name]]$orography, - var = 'orography', - longitude = values(list(lons.min, lons.max)), - longitude_reorder = circularsort, + orography_exp <- Start(dat = archive$System[[exp.name]]$orography, + var = 'orography', latitude = values(list(lats.min, lats.max)), latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = circularsort, transform = regrid_params$fcst.transform, - transform_params = list(grid = regrid_params$fcst.gridtype, - method = regrid_params$fcst.gridmethod), - transform_vars = c('latitude', 'longitude'), - return_vars = list(latitude = NULL, longitude = NULL), - synonims = list(longitude = c('lon', 'longitude'), - latitude = c('lat', 'latitude')), + transform_params = list(grid = regrid_params$fcst.gridtype, + method = regrid_params$fcst.gridmethod), + transform_vars = c('latitude','longitude'), + return_vars = list(lat = NULL, lon = NULL), + synonims = list(longitude = c('lon','longitude'), + latitude = c('lat','latitude')), num_procs = 1, retrieve = TRUE) - + if (attributes(orography_exp)$Variables$common$oro$unit != 'm') { + if (attributes(orography_exp)$Variables$common$oro$unit == 'km') { + orography_exp <- orography_exp * 1000 + attributes(orography_exp)$Variables$common$oro$unit <- 'm' + } else { + stop("Exp orography is not m or km, add option to orography_correction.R or fix the files") + } + } + + ## Set negative values to zero + orography_exp <- pmax(orography_exp, 0) + archive <- read_yaml("conf/archive_reference.yml")[[recipe$Run$filesystem]] ## load obs orography - orography_obs <- Start(dat = archive$Reference[[ref.name]]$orography, - var = 'orography', - longitude = values(list(lons.min, lons.max)), - longitude_reorder = circularsort, + orography_obs <- Start(dat = archive$Reference[[ref.name]]$orography, + var = 'orography', latitude = values(list(lats.min, lats.max)), latitude_reorder = Sort(), + longitude = values(list(lons.min, lons.max)), + longitude_reorder = circularsort, transform = regrid_params$obs.transform, transform_params = list(grid = regrid_params$obs.gridtype, method = regrid_params$obs.gridmethod), - transform_vars = c('latitude', 'longitude'), - return_vars = list(latitude = NULL, longitude = NULL), - synonims = list(longitude = c('lon', 'longitude'), - latitude = c('lat', 'latitude')), + transform_vars = c('latitude','longitude'), + return_vars = list(latitude = NULL, longitude = NULL), + synonims = list(longitude = c('lon','longitude'), + latitude = c('lat','latitude')), num_procs = 1, retrieve = TRUE) - - ## Set negative values to zero - orography_exp <- pmax(orography_exp, 0) - orography_obs <- pmax(orography_obs, 0) - - ## Calculate difference between orographies - oro_diff <- orography_exp - orography_obs - + if (attributes(orography_obs)$Variables$common$oro$unit != 'm') { + if (attributes(orography_obs)$Variables$common$oro$unit == 'km') { + orography_obs <- orography_obs * 1000 + attributes(orography_obs)$Variables$common$oro$unit <- 'm' + } else { + stop("Obs orography is not m or km, add option to orography_correction.R or fix the files") + } + } + # Convert to zero negative heights: + orography_obs[orography_obs < 0] <- 0 + orography_exp[orography_exp < 0] <- 0 + ## Apply lapse rate factor to correct temperature = -0.0065 K/m (-6.5 K/km) - oro_obs_corr <- oro_diff * -0.0065 - oro_obs_corr <- Reorder(data = drop(oro_obs_corr), order = c('latitude', 'longitude')) - - ## Apply correction to obs temperature data - data$obs$data <- Apply(data = list(data$obs$data, oro_obs_corr), - target_dims = c("latitude", "longitude"), - fun = "+", - ncores = recipe$Analysis$ncores)$output1 - data$obs$data <- Reorder(data$obs$data, names(data$obs$dims)) - + data$hcst$data <- Apply(list(data$hcst$data), + target_dims = c('latitude', 'longitude'), + fun = function(x, corr_const, + oro_obs, oro_exp) { + x + corr_const/1000 * + (oro_obs - oro_exp)}, + corr_const = correction_constant, + oro_obs = drop(orography_obs), + oro_exp = drop(orography_exp))$output1 + data$hcst$data <- Reorder(data$hcst$data, names(dim(data$obs$data))) + if (!is.null(data$fcst$data)) { + data$fcst$data <- Apply(list(data$fcst$data), + target_dims = c('latitude', 'longitude'), + fun = function(x, corr_const, + oro_obs, oro_exp) { + x + corr_const/1000 * + (oro_obs - oro_exp)}, + corr_const = correction_constant, + oro_obs = drop(orography_obs), + oro_exp = drop(orography_exp))$output1 + data$fcst$data <- Reorder(data$fcst$data, names(dim(data$obs$data))) + } return(data) + } -- GitLab From 6ca6a826ebb2373648b7b68981c59a5103f12aaf Mon Sep 17 00:00:00 2001 From: vagudets Date: Fri, 28 Mar 2025 12:28:06 +0100 Subject: [PATCH 2/2] Use get_archive() --- modules/Loading/R/orography_correction.R | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/modules/Loading/R/orography_correction.R b/modules/Loading/R/orography_correction.R index aa1ae035..190f32f9 100644 --- a/modules/Loading/R/orography_correction.R +++ b/modules/Loading/R/orography_correction.R @@ -1,4 +1,3 @@ -## TODO: remove paths to personal scratchs source("modules/Loading/R/get_regrid_params.R") source("modules/Loading/R/check_latlon.R") @@ -6,7 +5,7 @@ source("modules/Loading/R/check_latlon.R") # This correction is relevant for the mean_bias metric. # correction_constant is negative # -6.5ºK/km(ºC/km for the environmental lapse rate correction -# -9.8ºK/km (�C/km) for the adiabatic correction +# -9.8ºK/km (ºC/km) for the adiabatic correction orography_correction <- function(recipe, data, correction_constant = -6.5) { ref.name <- recipe$Analysis$Datasets$Reference$name @@ -16,7 +15,7 @@ orography_correction <- function(recipe, data, correction_constant = -6.5) { lons.min <- recipe$Analysis$Region$lonmin lons.max <- recipe$Analysis$Region$lonmax - archive <- read_yaml("conf/archive_seasonal.yml")[[recipe$Run$filesystem]] + archive <- get_archive(recipe) # Define regrid parameters: regrid_params <- get_regrid_params(recipe, archive) @@ -49,7 +48,6 @@ orography_correction <- function(recipe, data, correction_constant = -6.5) { ## Set negative values to zero orography_exp <- pmax(orography_exp, 0) - archive <- read_yaml("conf/archive_reference.yml")[[recipe$Run$filesystem]] ## load obs orography orography_obs <- Start(dat = archive$Reference[[ref.name]]$orography, var = 'orography', @@ -101,6 +99,5 @@ orography_correction <- function(recipe, data, correction_constant = -6.5) { data$fcst$data <- Reorder(data$fcst$data, names(dim(data$obs$data))) } return(data) - } -- GitLab