Transform difference when the region reaches the border
The transform results are slightly different between startR v1.0.0 and v0.1.4 when the selected region reaches one of the borders.
The transformation function in startR is CDORemapper, which uses s2dverification::CDORemap inside. And CDORemap uses cdo inside to do the interpolation. There are two possible ways to use cdo: (1) Use the global domain as the input, and use 'sellonlatbox' to crop the region (2) Use the region as the input (including the additional grid point assigned in 'transform_extra_cells'.
In the old version, when "latitude/longitude + transform_extra_cells" exceeds any of the two boundaries, it uses the whole domain to do the transformation then crops the selected region (method 1). In the new version, if "latitude/longitude + transform_extra_cells" exceeds the boundaries, it uses the indices until the boundary for transformation (method 2). These two methods cause a small difference at the border.
Take region lat = [0, 10]; lon = [0, 10]
as an example. lon = 0 is at the boundary when the region is defined as [0, 360] in CircularSort(). In the old version, the longitude index for cdo (in s2dverification::CDORemap, which is called by startR::CDORemapper) is the full domain, [1:1280]. In the new version, the longitude index for cdo is [1:38], 1 is the point lon = 0 and 38 is 36 (the grid point lon = 10) + 2 (transform_extra_cells = 2). The retrieved values are different at lon = 1.
Notice that:
- If the longitude is defined as [-180, 180], it doesn't reach the boundary so the results are identical between the old and new versions.
- When the region does not reach the boundary, both versions use "latitude/longitude + transform_extra_cells" for cdo (method 2). Therefore, the new version is more consistent regarding the transformation method.
- Load() returns the same result as the old version.
This issue should be kept in mind when using transformation in Start() or Load(). For now, there is no option to choose which method to be used. If you need to have the options or have query about this issue, please leave the comment below.
obs_path <- '/esarchive/recon/ecmwf/era5/monthly_mean/$var$_f1h/$var$_$sdate$.nc'
lons.min <- 0
lons.max <- 10
lats.min <- 0
lats.max <- 10
obs <- Start(dat = obs_path,
var = 'sfcWind',
sdate = '201811',
time = 'all',
latitude = values(list(lats.min, lats.max)),
latitude_reorder = Sort(decreasing = T),
longitude = values(list(lons.min, lons.max)),
longitude_reorder = CircularSort(0, 360),
synonims = list(longitude = c('lon', 'longitude'),
latitude = c('lat', 'latitude')),
transform = CDORemapper,
transform_extra_cells = 2,
transform_params = list(grid = 'r360x181',
method = 'conservative',
crop = c(lons.min, lons.max, lats.min, lats.max)),
transform_vars = c('latitude', 'longitude'),
return_vars = list(time = NULL, latitude = 'dat', longitude = 'dat'),
retrieve = T)
print(obs[1, 1, 1, 1, , ]
# Compare with Load()
pobs <- paste0('/esarchive/recon/ecmwf/era5/monthly_mean/',
'$VAR_NAME$_f1h/$VAR_NAME$_$YEAR$$MONTH$.nc')
obs_load <- Load(var = 'sfcWind',
obs = list(list(path = pobs)),
sdates = '20181101',
leadtimemin = 1,
leadtimemax = 1,
output = 'lonlat',
grid = 'r360x181',
method = 'conservative',
storefreq = 'monthly',
latmin = lats.min,
latmax = lats.max,
lonmin = lons.min,
lonmax = lons.max)
print(obs_load[1, 1, 1, 1, , ]