Skip to content
GitLab
Projects
Groups
Topics
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Sign in
Toggle navigation
Menu
Earth Sciences
CSDownscale
Compare revisions
62479ca38e4c5426681972661da0e1052821938e...72405abd8002fae8a8061b25bf98f43cdb320853
Commits (2)
adapted to hindcast-forecast
· 89a7953c
Jaume Ramon
authored
Nov 29, 2023
89a7953c
corrected coords by attrs for lat lon values
· 72405abd
Jaume Ramon
authored
Nov 29, 2023
72405abd
Hide whitespace changes
Inline
Side-by-side
R/Analogs.R
View file @
72405abd
...
...
@@ -99,9 +99,9 @@ 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
$
coord
s
[[
lat_dim
]],
exp_lons
=
exp
$
coord
s
[[
lon_dim
]],
obs_lats
=
obs
$
coord
s
[[
lat_dim
]],
obs_lons
=
obs
$
coord
s
[[
lon_dim
]],
grid_exp
=
grid_exp
,
nanalogs
=
nanalogs
,
res
<-
Analogs
(
exp
=
exp
$
data
,
obs
=
obs
$
data
,
exp_lats
=
exp
$
attr
s
[[
lat_dim
]],
exp_lons
=
exp
$
attr
s
[[
lon_dim
]],
obs_lats
=
obs
$
attr
s
[[
lat_dim
]],
obs_lons
=
obs
$
attr
s
[[
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
,
...
...
@@ -110,13 +110,13 @@ CST_Analogs <- function(exp, obs, grid_exp, obs2 = NULL, nanalogs = 3, fun_analo
# Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data
exp
$
data
<-
res
$
data
exp
$
dims
<-
dim
(
exp
$
data
)
exp
$
coord
s
[[
lon_dim
]]
<-
res
$
lon
exp
$
coord
s
[[
lat_dim
]]
<-
res
$
lat
exp
$
attr
s
[[
lon_dim
]]
<-
res
$
lon
exp
$
attr
s
[[
lat_dim
]]
<-
res
$
lat
obs
$
data
<-
res
$
obs
obs
$
dims
<-
dim
(
obs
$
data
)
obs
$
coord
s
[[
lon_dim
]]
<-
res
$
lon
obs
$
coord
s
[[
lat_dim
]]
<-
res
$
lat
obs
$
attr
s
[[
lon_dim
]]
<-
res
$
lon
obs
$
attr
s
[[
lat_dim
]]
<-
res
$
lat
res_s2dv
<-
list
(
exp
=
exp
,
obs
=
obs
)
return
(
res_s2dv
)
...
...
R/Intbc.R
View file @
72405abd
...
...
@@ -20,6 +20,10 @@
#'@param obs an 's2dv object' containing the observational field. The object
#'must have, at least, the dimensions latitude, longitude and start date. The object is
#'expected to be already subset for the desired region.
#'@param exp_cor an optional 's2dv_cube' object with named dimensions containing the seasonal
#'forecast experiment data. If the forecast is provided, it will be downscaled using the
#'hindcast and observations; if not provided, the hindcast will be downscaled instead. The
#'default value is NULL.
#'@param target_grid a character vector indicating the target grid to be passed to CDO.
#'It must be a grid recognised by CDO or a NetCDF file.
#'@param bc_method a character vector indicating the bias adjustment method to be applied after
...
...
@@ -67,7 +71,7 @@
#'res <- CST_Intbc(exp = exp, obs = obs, target_grid = 'r1280x640', bc_method = 'simple_bias', int_method = 'conservative')
#'@export
CST_Intbc
<-
function
(
exp
,
obs
,
target_grid
,
bc_method
,
int_method
=
NULL
,
points
=
NULL
,
CST_Intbc
<-
function
(
exp
,
obs
,
exp_cor
=
NULL
,
target_grid
,
bc_method
,
int_method
=
NULL
,
points
=
NULL
,
method_point_interp
=
NULL
,
lat_dim
=
"lat"
,
lon_dim
=
"lon"
,
sdate_dim
=
"sdate"
,
member_dim
=
"member"
,
region
=
NULL
,
ncores
=
NULL
,
...
)
{
...
...
@@ -79,25 +83,36 @@ 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
$
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
,
res
<-
Intbc
(
exp
=
exp
$
data
,
obs
=
obs
$
data
,
exp_cor
=
exp_cor
$
data
,
exp_lats
=
exp
$
attrs
[[
lat_dim
]],
exp_lons
=
exp
$
attrs
[[
lon_dim
]],
obs_lats
=
obs
$
attrs
[[
lat_dim
]],
obs_lons
=
obs
$
attrs
[[
lon_dim
]],
target_grid
=
target_grid
,
int_method
=
int_method
,
bc_method
=
bc_method
,
points
=
points
,
source_file_exp
=
exp
$
attrs
$
source_files
[
1
],
source_file_obs
=
obs
$
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
,
region
=
region
,
ncores
=
ncores
,
...
)
# Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data
exp
$
data
<-
res
$
data
exp
$
dims
<-
dim
(
exp
$
data
)
exp
$
coords
[[
lon_dim
]]
<-
res
$
lon
exp
$
coords
[[
lat_dim
]]
<-
res
$
lat
obs
$
data
<-
res
$
obs
obs
$
dims
<-
dim
(
obs
$
data
)
obs
$
coords
[[
lon_dim
]]
<-
res
$
lon
obs
$
coords
[[
lat_dim
]]
<-
res
$
lat
obs
$
attrs
[[
lon_dim
]]
<-
res
$
lon
obs
$
attrs
[[
lat_dim
]]
<-
res
$
lat
if
(
is.null
(
exp_cor
))
{
exp
$
data
<-
res
$
data
exp
$
dims
<-
dim
(
exp
$
data
)
exp
$
attrs
[[
lon_dim
]]
<-
res
$
lon
exp
$
attrs
[[
lat_dim
]]
<-
res
$
lat
res_s2dv
<-
list
(
exp
=
exp
,
obs
=
obs
)
}
else
{
exp_cor
$
data
<-
res
$
data
exp_cor
$
dims
<-
dim
(
exp_cor
$
data
)
exp_cor
$
attrs
[[
lon_dim
]]
<-
res
$
lon
exp_cor
$
attrs
[[
lat_dim
]]
<-
res
$
lat
res_s2dv
<-
list
(
exp
=
exp_cor
,
obs
=
obs
)
}
res_s2dv
<-
list
(
exp
=
exp
,
obs
=
obs
)
return
(
res_s2dv
)
}
...
...
@@ -123,7 +138,11 @@ CST_Intbc <- function(exp, obs, target_grid, bc_method, int_method = NULL, point
#''region'.
#'@param obs an array with named dimensions containing the observational field. The object
#'must have, at least, the dimensions latitude, longitude and start date. The object is
#'expected to be already subset for the desired region.
#'expected to be already subset for the desired region.
#'@param exp_cor an optional array with named dimensions containing the seasonal forecast
#'experiment data. If the forecast is provided, it will be downscaled using the hindcast and
#'observations; if not provided, the hindcast will be downscaled instead. The default value
#'is NULL.
#'@param exp_lats a numeric vector containing the latitude values in 'exp'. Latitudes must
#'range from -90 to 90.
#'@param exp_lons a numeric vector containing the longitude values in 'exp'. Longitudes
...
...
R/Interpolation.R
View file @
72405abd
...
...
@@ -68,7 +68,7 @@ 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
$
coord
s
[[
lat_dim
]],
lons
=
exp
$
coord
s
[[
lon_dim
]],
res
<-
Interpolation
(
exp
=
exp
$
data
,
lats
=
exp
$
attr
s
[[
lat_dim
]],
lons
=
exp
$
attr
s
[[
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
,
ncores
=
ncores
)
...
...
@@ -76,8 +76,8 @@ CST_Interpolation <- function(exp, points = NULL, method_remap = NULL, target_gr
# Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data
exp
$
data
<-
res
$
data
exp
$
dims
<-
dim
(
exp
$
data
)
exp
$
coord
s
[[
lon_dim
]]
<-
res
$
lon
exp
$
coord
s
[[
lat_dim
]]
<-
res
$
lat
exp
$
attr
s
[[
lon_dim
]]
<-
res
$
lon
exp
$
attr
s
[[
lat_dim
]]
<-
res
$
lat
res_s2dv
<-
list
(
exp
=
exp
,
obs
=
NULL
)
return
(
res_s2dv
)
...
...
R/Intlr.R
View file @
72405abd
...
...
@@ -22,6 +22,10 @@
#'@param obs an 's2dv object' containing the observational field. The object
#'must have, at least, the dimensions latitude, longitude and start date. The object is
#'expected to be already subset for the desired region.
#'@param exp_cor an optional 's2dv_cube' object with named dimensions containing the seasonal
#'forecast experiment data. If the forecast is provided, it will be downscaled using the
#'hindcast and observations; if not provided, the hindcast will be downscaled instead. The
#'default value is NULL.
#'@param lr_method a character vector indicating the linear regression method to be applied
#'after the interpolation. Accepted methods are 'basic', 'large-scale' and '4nn'. The 'basic'
#'method fits a linear regression using high resolution observations as predictands and the
...
...
@@ -86,7 +90,7 @@
#'obs <- s2dv_cube(data = obs, lat = obs_lats, lon = obs_lons)
#'res <- CST_Intlr(exp = exp, obs = obs, target_grid = 'r1280x640', lr_method = 'basic', int_method = 'conservative')
#'@export
CST_Intlr
<-
function
(
exp
,
obs
,
lr_method
,
target_grid
=
NULL
,
points
=
NULL
,
int_method
=
NULL
,
CST_Intlr
<-
function
(
exp
,
obs
,
exp_cor
=
NULL
,
lr_method
,
target_grid
=
NULL
,
points
=
NULL
,
int_method
=
NULL
,
method_point_interp
=
NULL
,
predictors
=
NULL
,
lat_dim
=
"lat"
,
lon_dim
=
"lon"
,
sdate_dim
=
"sdate"
,
time_dim
=
"time"
,
member_dim
=
"member"
,
large_scale_predictor_dimname
=
'vars'
,
loocv
=
TRUE
,
region
=
NULL
,
ncores
=
NULL
)
{
...
...
@@ -98,28 +102,54 @@ CST_Intlr <- function(exp, obs, lr_method, target_grid = NULL, points = NULL, in
if
(
!
inherits
(
obs
,
's2dv_cube'
))
{
stop
(
"Parameter 'obs' must be of the class 's2dv_cube'"
)
}
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
],
if
(
identical
(
lr_method
,
'basic'
)
|
identical
(
lr_method
,
'4nn'
))
{
if
(
!
inherits
(
exp_cor
,
's2dv_cube'
))
{
stop
(
"Parameter 'exp_cor' must be of the class 's2dv_cube'"
)
}
exp_cor_aux
<-
exp_cor
$
data
# when large-scale is selected, the forecast object does not have to be an s2dv_cube object
}
else
if
(
identical
(
lr_method
,
'large-scale'
))
{
exp_cor_aux
<-
exp_cor
}
res
<-
Intlr
(
exp
=
exp
$
data
,
obs
=
obs
$
data
,
exp_cor
=
exp_cor_aux
,
exp_lats
=
exp
$
attrs
[[
lat_dim
]],
exp_lons
=
exp
$
attrs
[[
lon_dim
]],
obs_lats
=
obs
$
attrs
[[
lat_dim
]],
obs_lons
=
obs
$
attrs
[[
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
,
lat_dim
=
lat_dim
,
lon_dim
=
lon_dim
,
sdate_dim
=
sdate_dim
,
time_dim
=
time_dim
,
member_dim
=
member_dim
,
large_scale_predictor_dimname
=
large_scale_predictor_dimname
,
loocv
=
loocv
,
region
=
region
,
ncores
=
ncores
)
# Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data
exp
$
data
<-
res
$
data
exp
$
dims
<-
dim
(
exp
$
data
)
exp
$
coords
[[
lon_dim
]]
<-
res
$
lon
exp
$
coords
[[
lat_dim
]]
<-
res
$
lat
obs
$
data
<-
res
$
obs
obs
$
dims
<-
dim
(
obs
$
data
)
obs
$
coords
[[
lon_dim
]]
<-
res
$
lon
obs
$
coords
[[
lat_dim
]]
<-
res
$
lat
res_s2dv
<-
list
(
exp
=
exp
,
obs
=
obs
)
obs
$
attrs
[[
lon_dim
]]
<-
res
$
lon
obs
$
attrs
[[
lat_dim
]]
<-
res
$
lat
if
(
is.null
(
exp_cor
))
{
exp
$
data
<-
res
$
data
exp
$
dims
<-
dim
(
exp
$
data
)
exp
$
attrs
[[
lon_dim
]]
<-
res
$
lon
exp
$
attrs
[[
lat_dim
]]
<-
res
$
lat
res_s2dv
<-
list
(
exp
=
exp
,
obs
=
obs
)
}
else
{
if
(
identical
(
lr_method
,
'basic'
)
|
identical
(
lr_method
,
'4nn'
))
{
exp_cor
$
data
<-
res
$
data
exp_cor
$
dims
<-
dim
(
exp_cor
$
data
)
exp_cor
$
attrs
[[
lon_dim
]]
<-
res
$
lon
exp_cor
$
attrs
[[
lat_dim
]]
<-
res
$
lat
# when large-scale is selected, the forecast object does not have to be an s2dv_cube object
}
else
if
(
identical
(
lr_method
,
'large-scale'
))
{
exp_cor
<-
suppressWarnings
(
s2dv_cube
(
res
$
data
,
lat
=
res
$
lat
,
lon
=
res
$
lon
))
}
res_s2dv
<-
list
(
exp
=
exp_cor
,
obs
=
obs
)
}
return
(
res_s2dv
)
}
...
...
@@ -147,6 +177,10 @@ CST_Intlr <- function(exp, obs, lr_method, target_grid = NULL, points = NULL, in
#'@param obs an array with named dimensions containing the observational field. The object
#'must have, at least, the dimensions latitude, longitude and start date. The object is
#'expected to be already subset for the desired region.
#'@param exp_cor an optional array with named dimensions containing the seasonal forecast
#'experiment data. If the forecast is provided, it will be downscaled using the hindcast and
#'observations; if not provided, the hindcast will be downscaled instead. The default value
#'is NULL.
#'@param exp_lats a numeric vector containing the latitude values in 'exp'. Latitudes must
#'range from -90 to 90.
#'@param exp_lons a numeric vector containing the longitude values in 'exp'. Longitudes
...
...
@@ -222,7 +256,7 @@ CST_Intlr <- function(exp, obs, lr_method, target_grid = NULL, points = NULL, in
#'res <- Intlr(exp = exp, obs = obs, exp_lats = exp_lats, exp_lons = exp_lons, obs_lats = obs_lats,
#'obs_lons = obs_lons, target_grid = 'r1280x640', lr_method = 'basic', int_method = 'conservative')
#'@export
Intlr
<-
function
(
exp
,
obs
,
exp_lats
,
exp_lons
,
obs_lats
,
obs_lons
,
lr_method
,
target_grid
=
NULL
,
points
=
NULL
,
Intlr
<-
function
(
exp
,
obs
,
exp_cor
=
NULL
,
exp_lats
,
exp_lons
,
obs_lats
,
obs_lons
,
lr_method
,
target_grid
=
NULL
,
points
=
NULL
,
int_method
=
NULL
,
method_point_interp
=
NULL
,
source_file_exp
=
NULL
,
source_file_obs
=
NULL
,
predictors
=
NULL
,
lat_dim
=
"lat"
,
lon_dim
=
"lon"
,
sdate_dim
=
"sdate"
,
time_dim
=
"time"
,
member_dim
=
"member"
,
region
=
NULL
,
large_scale_predictor_dimname
=
'vars'
,
loocv
=
TRUE
,
...
...
@@ -274,6 +308,35 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t
"'sdate_dim'"
)
}
if
(
!
is.null
(
exp_cor
))
{
if
(
is.na
(
match
(
sdate_dim
,
names
(
dim
(
exp_cor
)))))
{
stop
(
"Missing start date dimension in 'exp_cor', or does not match the parameter "
,
"'sdate_dim'"
)
}
if
(
is.na
(
match
(
member_dim
,
names
(
dim
(
exp_cor
)))))
{
stop
(
"Missing member dimension in 'exp_cor', or does not match the parameter 'member_dim'"
)
}
if
(
loocv
)
{
# loocv equal to false to train with the whole hindcast and predict with the forecast
loocv
<-
FALSE
}
if
(
is.null
(
predictors
))
{
if
(
is.na
(
match
(
lon_dim
,
names
(
dim
(
exp_cor
)))))
{
stop
(
"Missing longitude dimension in 'exp_cor', or does not match the parameter "
,
"'lon_dim'"
)
}
if
(
is.na
(
match
(
lat_dim
,
names
(
dim
(
exp_cor
)))))
{
stop
(
"Missing latitude dimension in 'exp_cor', or does not match the parameter "
,
"'lat_dim'"
)
}
}
}
# When observations are pointwise
if
(
!
is.null
(
points
)
&
!
is.na
(
match
(
"location"
,
names
(
dim
(
obs
)))))
{
point_obs
<-
T
...
...
@@ -295,10 +358,6 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t
"through the parameter 'method_point_interp'."
)
}
# sdate must be the time dimension in the input data
stopifnot
(
sdate_dim
%in%
names
(
dim
(
exp
)))
stopifnot
(
sdate_dim
%in%
names
(
dim
(
obs
)))
# the code is not yet prepared to handle members in the observations
restore_ens
<-
FALSE
if
(
member_dim
%in%
names
(
dim
(
obs
)))
{
...
...
@@ -321,6 +380,17 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t
stopifnot
(
sdate_dim
%in%
names
(
dim
(
predictors
)))
stopifnot
(
dim
(
predictors
)[
sdate_dim
]
==
dim
(
exp
)[
sdate_dim
])
}
# forecasts for the large scale predictors
if
(
!
is.null
(
exp_cor
))
{
if
(
is.na
(
match
(
large_scale_predictor_dimname
,
names
(
dim
(
exp_cor
)))))
{
stop
(
"Missing large scale predictor dimension in 'exp_cor', or does not match the parameter "
,
"'large_scale_predictor_dimname'"
)
}
if
(
!
identical
(
dim
(
exp_cor
)[
names
(
dim
(
exp_cor
))
==
large_scale_predictor_dimname
],
dim
(
predictors
)[
names
(
dim
(
predictors
))
==
large_scale_predictor_dimname
]))
{
stop
(
"Large scale predictor dimension in exp_cor and predictors must be identical."
)
}
}
}
## ncores
if
(
!
is.null
(
ncores
))
{
...
...
@@ -351,6 +421,13 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t
points
=
points
,
method_point_interp
=
method_point_interp
,
source_file
=
source_file_exp
,
lat_dim
=
lat_dim
,
lon_dim
=
lon_dim
,
method_remap
=
int_method
,
region
=
region
,
ncores
=
ncores
)
if
(
!
is.null
(
exp_cor
)
&
is.null
(
predictors
))
{
exp_cor_interpolated
<-
Interpolation
(
exp
=
exp_cor
,
lats
=
exp_lats
,
lons
=
exp_lons
,
target_grid
=
target_grid
,
points
=
points
,
method_point_interp
=
method_point_interp
,
source_file
=
source_file_exp
,
lat_dim
=
lat_dim
,
lon_dim
=
lon_dim
,
method_remap
=
int_method
,
region
=
region
,
ncores
=
ncores
)
}
# If after interpolating 'exp' data the coordinates do not match, the obs data is interpolated to
# the same grid to force the matching
...
...
@@ -383,6 +460,17 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t
target_dims_predictor
<-
sdate_dim
target_dims_predictand
<-
sdate_dim
if
(
!
is.null
(
exp_cor
))
{
aux_dim
<-
NULL
forecast
<-
exp_cor_interpolated
$
data
target_dims_predictor
<-
c
(
sdate_dim
,
member_dim
)
target_dims_forecast
<-
c
(
sdate_dim
,
member_dim
)
}
else
{
forecast
<-
NULL
target_dims_forecast
<-
NULL
}
}
# (Multi) linear regression with large-scale predictors
...
...
@@ -398,6 +486,13 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t
target_dims_predictor
<-
c
(
sdate_dim
,
large_scale_predictor_dimname
)
target_dims_predictand
<-
sdate_dim
if
(
!
is.null
(
exp_cor
))
{
aux_dim
<-
large_scale_predictor_dimname
forecast
<-
exp_cor
target_dims_predictor
<-
c
(
sdate_dim
,
large_scale_predictor_dimname
,
member_dim
)
target_dims_forecast
<-
c
(
sdate_dim
,
large_scale_predictor_dimname
,
member_dim
)
}
}
# Multi-linear regression with the four nearest neighbours
...
...
@@ -406,8 +501,14 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t
else
if
(
lr_method
==
'4nn'
)
{
predictor
<-
.find_nn
(
coar
=
exp
,
lats_hres
=
obs_lats
,
lons_hres
=
obs_lons
,
lats_coar
=
exp_lats
,
lons_coar
=
exp_lons
,
lat_dim
=
lat_dim
,
lon_dim
=
lon_dim
,
nn
=
4
,
ncores
=
ncores
)
lons_coar
=
exp_lons
,
lat_dim
=
lat_dim
,
lon_dim
=
lon_dim
,
nn
=
4
,
ncores
=
ncores
)
if
(
!
is.null
(
exp_cor
))
{
aux_dim
<-
'nn'
forecast
<-
.find_nn
(
coar
=
exp_cor
,
lats_hres
=
obs_lats
,
lons_hres
=
obs_lons
,
lats_coar
=
exp_lats
,
lons_coar
=
exp_lons
,
lat_dim
=
lat_dim
,
lon_dim
=
lon_dim
,
nn
=
4
,
ncores
=
ncores
)
}
if
(
is.null
(
points
)
|
(
"location"
%in%
names
(
dim
(
obs
))))
{
if
(
!
is.null
(
target_grid
))
{
warning
(
"Interpolating to the 'obs' grid"
)
...
...
@@ -433,27 +534,34 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t
predictor
<-
predictor
$
data
predictand
<-
predictand
$
data
}
target_dims_predictor
<-
c
(
sdate_dim
,
'nn'
)
target_dims_predictand
<-
sdate_dim
target_dims_predictand
<-
sdate_dim
if
(
!
is.null
(
exp_cor
))
{
target_dims_predictor
<-
c
(
sdate_dim
,
'nn'
,
member_dim
)
target_dims_forecast
<-
c
(
sdate_dim
,
'nn'
,
member_dim
)
}
else
{
target_dims_predictor
<-
c
(
sdate_dim
,
'nn'
)
}
}
else
{
stop
(
paste0
(
lr_method
,
" method is not implemented yet"
))
}
print
(
paste0
(
'dim predictor'
,
dim
(
predictor
)))
print
(
paste0
(
'dim predictand'
,
dim
(
predictand
)))
print
(
dim
(
list
(
predictor
[
1
])))
# Apply the linear regressions
res
<-
Apply
(
list
(
predictor
,
predictand
),
target_dims
=
list
(
target_dims_predictor
,
target_dims_predictand
),
fun
=
.intlr
,
loocv
=
loocv
,
ncores
=
ncores
)
$
output1
names
(
dim
(
res
))[
1
]
<-
sdate_dim
# names(dim(res))[which(names(dim(res)) == '')]
# Apply the linear regressions
# case hindcast - forecast
if
(
!
is.null
(
exp_cor
))
{
res
<-
Apply
(
list
(
predictor
,
predictand
,
forecast
),
target_dims
=
list
(
target_dims_predictor
,
target_dims_predictand
,
target_dims_forecast
),
fun
=
.intlr
,
loocv
=
loocv
,
aux_dim
=
aux_dim
,
ncores
=
ncores
)
$
output1
}
# case hindcast only
else
{
res
<-
Apply
(
list
(
predictor
,
predictand
),
target_dims
=
list
(
target_dims_predictor
,
target_dims_predictand
),
fun
=
.intlr
,
loocv
=
loocv
,
ncores
=
ncores
)
$
output1
names
(
dim
(
res
))[
1
]
<-
sdate_dim
}
# restore ensemble dimension in observations if it existed originally
if
(
restore_ens
)
{
...
...
@@ -463,30 +571,56 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t
# Return a list of three elements
res
<-
list
(
data
=
res
,
obs
=
predictand
,
lon
=
lons
,
lat
=
lats
)
# for testing 4nn
#x <- predictor[10,10,1,1,1,,,]
#x <- Reorder(x,c(2,3,1))
#y <- predictand[1,1,1,10,,10]
#f <- forecast[10,10,1,1,1,,,]
#f <- Reorder(f,c(2,1))
#f <- InsertDim(f,1,1,'sdate')
# large-scale
#x <- Reorder(predictor,c(1,3,2))
#y <- predictand[1,1,1,10,,10]
#f <- Reorder(forecast,c(1,3,2))
return
(
res
)
}
#-----------------------------------
# Atomic function to generate and apply the linear regressions
#-----------------------------------
.intlr
<-
function
(
x
,
y
,
loocv
)
{
tmp_df
<-
data.frame
(
x
=
x
,
y
=
y
)
.intlr
<-
function
(
x
,
y
,
f
=
NULL
,
loocv
,
aux_dim
=
NULL
)
{
if
(
!
is.null
(
f
))
{
if
(
!
is.null
(
aux_dim
))
{
tmp_df
<-
data.frame
(
x
=
adply
(
x
,
.margins
=
3
,
.id
=
NULL
),
y
=
y
)
}
else
{
tmp_df
<-
data.frame
(
x
=
as.vector
(
x
),
y
=
y
)
}
}
else
{
tmp_df
<-
data.frame
(
x
=
x
,
y
=
y
)
}
# if the data is all NA, force return return NA
if
(
all
(
is.na
(
tmp_df
))
|
(
sum
(
apply
(
tmp_df
,
2
,
function
(
x
)
!
all
(
is.na
(
x
))))
==
1
))
{
n
<-
nrow
(
tmp_df
)
res
<-
rep
(
NA
,
n
)
if
(
is.null
(
f
))
{
res
<-
rep
(
NA
,
nrow
(
tmp_df
))
}
else
{
if
(
!
is.null
(
aux_dim
))
{
res
<-
array
(
NA
,
dim
(
f
)[
names
(
dim
(
f
))
!=
aux_dim
])
}
else
{
res
<-
array
(
NA
,
dim
(
f
))
}
}
}
else
{
# training
lm1
<-
.train_lm
(
df
=
tmp_df
,
loocv
=
loocv
)
# prediction
res
<-
.pred_lm
(
lm1
=
lm1
,
df
=
tmp_df
,
loocv
=
loocv
)
res
<-
.pred_lm
(
lm1
=
lm1
,
df
=
tmp_df
,
f
=
f
,
loocv
=
loocv
,
aux_dim
=
aux_dim
)
}
return
(
res
)
return
(
res
)
}
#-----------------------------------
...
...
@@ -517,7 +651,7 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t
#-----------------------------------
# Function to apply the linear regressions.
#-----------------------------------
.pred_lm
<-
function
(
df
,
lm1
,
loocv
)
{
.pred_lm
<-
function
(
df
,
lm1
,
f
,
loocv
,
aux_dim
)
{
if
(
loocv
)
{
pred_vals
<-
sapply
(
1
:
nrow
(
df
),
function
(
j
)
{
...
...
@@ -527,13 +661,46 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t
return
(
predict
(
lm1
[[
j
]],
df
[
j
,]))
}})
}
else
{
if
(
!
is.na
(
lm1
))
{
pred_vals_ls
<-
lapply
(
lm1
,
predict
,
data
=
df
)
pred_vals
<-
unlist
(
pred_vals_ls
)
if
(
!
is.na
(
lm1
))
{
# case to downscale hindcasts
if
(
is.null
(
f
))
{
pred_vals_ls
<-
lapply
(
lm1
,
predict
,
data
=
df
)
pred_vals
<-
unlist
(
pred_vals_ls
)
# case to downscale forecasts
}
else
{
if
(
!
is.null
(
aux_dim
))
{
fcst_df
<-
adply
(
f
,
.margins
=
3
,
.id
=
NULL
)
# 4nn
if
(
identical
(
aux_dim
,
'nn'
))
{
colnames
(
fcst_df
)
<-
paste0
(
'x.'
,
1
:
4
)
pred_vals_ls
<-
lapply
(
lm1
,
predict
,
newdata
=
fcst_df
)
pred_vals
<-
unlist
(
pred_vals_ls
)
dim
(
pred_vals
)
<-
dim
(
f
)[
names
(
dim
(
f
))
!=
aux_dim
]
}
else
{
pred_vals_ls
<-
lapply
(
lm1
,
predict
,
newdata
=
data.frame
(
x
=
fcst_df
))
pred_vals
<-
unlist
(
pred_vals_ls
)
dim
(
pred_vals
)
<-
dim
(
f
)[
names
(
dim
(
f
))
!=
aux_dim
]
}
# basic
}
else
{
pred_vals_ls
<-
lapply
(
lm1
,
predict
,
newdata
=
data.frame
(
x
=
as.vector
(
f
)))
pred_vals
<-
unlist
(
pred_vals_ls
)
dim
(
pred_vals
)
<-
dim
(
f
)
}
}
}
else
{
pred_vals
<-
rep
(
NA
,
nrow
(
df
))
if
(
is.null
(
f
))
{
pred_vals
<-
rep
(
NA
,
nrow
(
df
))
}
else
{
if
(
!
is.null
(
aux_dim
))
{
pred_vals
<-
array
(
NA
,
dim
(
f
)[
names
(
dim
(
f
))
!=
aux_dim
])
}
else
{
pred_vals
<-
array
(
NA
,
dim
(
f
))
}
}
}
}
}
return
(
pred_vals
)
}
...
...
R/LogisticReg.R
View file @
72405abd
...
...
@@ -99,9 +99,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
$
coord
s
[[
lat_dim
]],
exp_lons
=
exp
$
coord
s
[[
lon_dim
]],
obs_lats
=
obs
$
coord
s
[[
lat_dim
]],
obs_lons
=
obs
$
coord
s
[[
lon_dim
]],
target_grid
=
target_grid
,
res
<-
LogisticReg
(
exp
=
exp
$
data
,
obs
=
obs
$
data
,
exp_lats
=
exp
$
attr
s
[[
lat_dim
]],
exp_lons
=
exp
$
attr
s
[[
lon_dim
]],
obs_lats
=
obs
$
attr
s
[[
lat_dim
]],
obs_lons
=
obs
$
attr
s
[[
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,13 +112,13 @@ 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
$
dims
<-
dim
(
exp
$
data
)
exp
$
coord
s
[[
lon_dim
]]
<-
res
$
lon
exp
$
coord
s
[[
lat_dim
]]
<-
res
$
lat
exp
$
attr
s
[[
lon_dim
]]
<-
res
$
lon
exp
$
attr
s
[[
lat_dim
]]
<-
res
$
lat
obs
$
data
<-
res
$
obs
obs
$
dims
<-
dim
(
obs
$
data
)
obs
$
coord
s
[[
lon_dim
]]
<-
res
$
lon
obs
$
coord
s
[[
lat_dim
]]
<-
res
$
lat
obs
$
attr
s
[[
lon_dim
]]
<-
res
$
lon
obs
$
attr
s
[[
lat_dim
]]
<-
res
$
lat
res_s2dv
<-
list
(
exp
=
exp
,
obs
=
obs
)
return
(
res_s2dv
)
...
...