diff --git a/R/Analogs.R b/R/Analogs.R
index 274a4e095135325486549697cf75c3ff4a4e97a5..a69a66dfc80b0943bffcfc3085da5d25145edcff 100644
--- a/R/Analogs.R
+++ b/R/Analogs.R
@@ -99,21 +99,24 @@ CST_Analogs <- function(exp, obs, grid_exp, obs2 = NULL, nanalogs = 3, fun_analo
     stop("Parameter 'obs' must be of the class 's2dv_cube'")
   }
 
-  res <- Analogs(exp = exp$data, obs = obs$data, exp_lats = exp$attrs$lat, exp_lons = exp$attrs$lon,
-                 obs_lats = obs$attrs$lat, obs_lons = obs$attrs$lon, grid_exp = grid_exp, 
-                 nanalogs = nanalogs, fun_analog = fun_analog, lat_dim = lat_dim, lon_dim = lon_dim, 
+  res <- Analogs(exp = exp$data, obs = obs$data, exp_lats = exp$coords[[lat_dim]], 
+                 exp_lons = exp$coords[[lon_dim]], obs_lats = obs$coords[[lat_dim]], 
+                 obs_lons = obs$coords[[lon_dim]], grid_exp = grid_exp, nanalogs = nanalogs, 
+                 fun_analog = fun_analog, lat_dim = lat_dim, lon_dim = lon_dim, 
                  sdate_dim = sdate_dim, time_dim = time_dim, member_dim = member_dim, 
                  region = region, return_indices = return_indices, loocv_window = loocv_window, 
                  ncores = ncores)
 
   # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data
   exp$data <- res$data
-  exp$lon <- res$lon
-  exp$lat <- res$lat
+  exp$dims <- dim(exp$data)
+  exp$coords[[lon_dim]] <- res$lon
+  exp$coords[[lat_dim]] <- res$lat
 
   obs$data <- res$obs
-  obs$lat <- res$lat
-  obs$lon <- res$lon
+  obs$dims <- dim(obs$data)
+  obs$coords[[lon_dim]] <- res$lon
+  obs$coords[[lat_dim]] <- res$lat
 
   res_s2dv <- list(exp = exp, obs = obs)
   return(res_s2dv)
diff --git a/R/Intbc.R b/R/Intbc.R
index 20f892a86931f30584e5778f09a6da72601753ba..27de3735dbe88f0a60f8892406e287b55f40a0f7 100644
--- a/R/Intbc.R
+++ b/R/Intbc.R
@@ -78,8 +78,8 @@ CST_Intbc <- function(exp, obs, target_grid, bc_method, int_method = NULL, point
     stop("Parameter 'obs' must be of the class 's2dv_cube'")
   }
 
-  res <- Intbc(exp = exp$data, obs = obs$data, exp_lats = exp$attrs$lat, exp_lons = exp$attrs$lon, 
-               obs_lats = obs$attrs$lat, obs_lons = obs$attrs$lon, target_grid = target_grid, 
+  res <- Intbc(exp = exp$data, obs = obs$data, exp_lats = exp$coords[[lat_dim]], exp_lons = exp$coords[[lon_dim]], 
+               obs_lats = obs$coords[[lat_dim]], obs_lons = obs$coords[[lon_dim]], target_grid = target_grid, 
                int_method = int_method, bc_method = bc_method, points = points, 
                source_file = exp$attrs$source_files[1], method_point_interp = method_point_interp, 
                lat_dim = lat_dim, lon_dim = lon_dim, sdate_dim = sdate_dim, member_dim = member_dim, 
@@ -87,12 +87,14 @@ CST_Intbc <- function(exp, obs, target_grid, bc_method, int_method = NULL, point
 
   # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data
   exp$data <- res$data
-  exp$lon <- res$lon
-  exp$lat <- res$lat
-  
+  exp$dims <- dim(exp$data)
+  exp$coords[[lon_dim]] <- res$lon
+  exp$coords[[lat_dim]] <- res$lat
+
   obs$data <- res$obs
-  obs$lat <- res$lat
-  obs$lon <- res$lon
+  obs$dims <- dim(obs$data)
+  obs$coords[[lon_dim]] <- res$lon
+  obs$coords[[lat_dim]] <- res$lat
 
   res_s2dv <- list(exp = exp, obs = obs)
   return(res_s2dv)
diff --git a/R/Interpolation.R b/R/Interpolation.R
index 62c55949a98307cef65094e97e85fc4c60a09356..1599bf3be04dae19b40ef65466b03bfd913e1d2f 100644
--- a/R/Interpolation.R
+++ b/R/Interpolation.R
@@ -67,15 +67,16 @@ CST_Interpolation <- function(exp, points = NULL, method_remap = NULL, target_gr
     stop("The name of the latitude/longitude dimensions in 'exp$data' must match the parametres 'lat_dim' and 'lon_dim'")
   }
 
-  res <- Interpolation(exp = exp$data, lats = exp$attrs$lat, lons = exp$attrs$lon, 
+  res <- Interpolation(exp = exp$data, lats = exp$coords[[lat_dim]], lons = exp$coords[[lon_dim]], 
                        source_file = exp$attrs$source_files[1], points = points,
                        method_remap = method_remap, target_grid = target_grid, lat_dim = lat_dim,
                        lon_dim = lon_dim, region = region, method_point_interp = method_point_interp)
 
   # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data
   exp$data <- res$data
-  exp$lon <- res$lon
-  exp$lat <- res$lat
+  exp$dims <- dim(exp$data)
+  exp$coords[[lon_dim]] <- res$lon
+  exp$coords[[lat_dim]] <- res$lat
 
   res_s2dv <- list(exp = exp, obs  = NULL)
   return(res_s2dv)
diff --git a/R/Intlr.R b/R/Intlr.R
index 5de4ff4893aab7e1a6736cdca5b063384559e532..24c909f367b21b03e10b6c0e79a42a9483654f09 100644
--- a/R/Intlr.R
+++ b/R/Intlr.R
@@ -99,8 +99,8 @@ CST_Intlr <- function(exp, obs, lr_method, target_grid = NULL, points = NULL, in
     stop("Parameter 'obs' must be of the class 's2dv_cube'")
   }
 
-  res <- Intlr(exp = exp$data, obs = obs$data, exp_lats = exp$attrs$lat, exp_lons = exp$attrs$lon, 
-               obs_lats = obs$attrs$lat, obs_lons = obs$attrs$lon, points = points, 
+  res <- Intlr(exp = exp$data, obs = obs$data, exp_lats = exp$coords[[lat_dim]], exp_lons = exp$coords[[lon_dim]],
+               obs_lats = obs$coords[[lat_dim]], obs_lons = obs$coords[[lon_dim]], points = points, 
                source_file_exp = exp$attrs$source_files[1], source_file_obs = obs$attrs$source_files[1], 
                target_grid = target_grid, lr_method = lr_method, int_method = int_method, 
                method_point_interp = method_point_interp, predictors = predictors, 
@@ -110,12 +110,14 @@ CST_Intlr <- function(exp, obs, lr_method, target_grid = NULL, points = NULL, in
 
   # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data
   exp$data <- res$data
-  exp$lon <- res$lon
-  exp$lat <- res$lat
+  exp$dims <- dim(exp$data)
+  exp$coords[[lon_dim]] <- res$lon
+  exp$coords[[lat_dim]] <- res$lat
 
   obs$data <- res$obs
-  obs$lat <- res$lat
-  obs$lon <- res$lon
+  obs$dims <- dim(obs$data)
+  obs$coords[[lon_dim]] <- res$lon
+  obs$coords[[lat_dim]] <- res$lat
 
   res_s2dv <- list(exp = exp, obs = obs)
   return(res_s2dv)
diff --git a/R/LogisticReg.R b/R/LogisticReg.R
index 1f19712af153d7f3d853fced6fe5e82dbfa7822b..c514d2541c1c74c5e53d50b11ea9c1e4775ab6b5 100644
--- a/R/LogisticReg.R
+++ b/R/LogisticReg.R
@@ -101,8 +101,9 @@ CST_LogisticReg <- function(exp, obs, target_grid, int_method = NULL, log_reg_me
     stop("Parameter 'obs' must be of the class 's2dv_cube'")
   }
 
-  res <- LogisticReg(exp = exp$data, obs = obs$data, exp_lats = exp$attrs$lat, exp_lons = exp$attrs$lon, 
-                     obs_lats = obs$attrs$lat, obs_lons = obs$attrs$lon, target_grid = target_grid,
+  res <- LogisticReg(exp = exp$data, obs = obs$data, exp_lats = exp$coords[[lat_dim]], 
+                     exp_lons = exp$coords[[lon_dim]], obs_lats = obs$coords[[lat_dim]], 
+                     obs_lons = obs$coords[[lon_dim]], target_grid = target_grid,
                      probs_cat = probs_cat, return_most_likely_cat = return_most_likely_cat,
                      int_method = int_method, log_reg_method = log_reg_method, points = points, 
                      method_point_interp = method_point_interp, lat_dim = lat_dim, 
@@ -112,12 +113,14 @@ CST_LogisticReg <- function(exp, obs, target_grid, int_method = NULL, log_reg_me
 
   # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data
   exp$data <- res$data
-  exp$lon <- res$lon
-  exp$lat <- res$lat
+  exp$dims <- dim(exp$data)
+  exp$coords[[lon_dim]] <- res$lon
+  exp$coords[[lat_dim]] <- res$lat
 
   obs$data <- res$obs
-  obs$lat <- res$lat
-  obs$lon <- res$lon
+  obs$dims <- dim(obs$data)
+  obs$coords[[lon_dim]] <- res$lon
+  obs$coords[[lat_dim]] <- res$lat
 
   res_s2dv <- list(exp = exp, obs = obs)
   return(res_s2dv)