From 63e47e7db69088a8517148fec8665a40fea3bb7d Mon Sep 17 00:00:00 2001 From: aho Date: Tue, 2 Nov 2021 14:17:44 +0100 Subject: [PATCH] Update vignettes with s2dv --- vignettes/anomaly_agreement.Rmd | 30 ++++++++++++----------- vignettes/diurnaltemp.Rmd | 42 ++++++++++++++++++++------------- vignettes/extreme_indices.Rmd | 41 ++++++++++++++++---------------- vignettes/heatcoldwaves.Rmd | 16 +++++++------ 4 files changed, 71 insertions(+), 58 deletions(-) diff --git a/vignettes/anomaly_agreement.Rmd b/vignettes/anomaly_agreement.Rmd index 85ff031..a432a25 100644 --- a/vignettes/anomaly_agreement.Rmd +++ b/vignettes/anomaly_agreement.Rmd @@ -32,7 +32,7 @@ All the other R packages involved can be installed directly from CRAN and loaded ```r library(abind) -library(s2dverification) +library(s2dv) library(ggplot2) ``` @@ -170,18 +170,18 @@ model time lon lat 4 90 12 8 ``` -+ Now, the mean climatology value for each grid point and model can be computed by using the `Mean1Dim()` function belonging to the **s2dverification package**. The position of the temporal dimension should be specified in parameter `posdim`. ++ Now, the mean climatology value for each grid point and model can be computed by using the `MeanDims()` function belonging to the **s2dv package**. The position of the temporal dimension should be specified in parameter `posdim`. ```r -climatology <- Mean1Dim(summer_historical$data, posdim = 2) +climatology <- MeanDims(summer_historical$data, dims = 'time') ``` -+ `Season()` function from **s2dverification package** returns the mean annual time series for the selected season by defining the parameters of the initial month of the data (`monini = 1`), the first month of the season (`moninf = 6`) and the final month of the season (`monsup = 8`). The last two parameters, `moninf` and `monsup`, have their origin with respect to the first one `monini`. ++ `Season()` function from **s2dv package** returns the mean annual time series for the selected season by defining the parameters of the initial month of the data (`monini = 1`), the first month of the season (`moninf = 6`) and the final month of the season (`monsup = 8`). The last two parameters, `moninf` and `monsup`, have their origin with respect to the first one `monini`. ```r -summer_projection <- Season(multimodel_projection, posdim = 3, +summer_projection <- Season(multimodel_projection, time_dim = 'time', monini = 1, moninf = 6, monsup = 8) ``` @@ -192,15 +192,17 @@ By running the next lines, it is possible to check the dimensions of the data: model lon lat 4 12 8 > dim(summer_projection) -model var time lon lat - 4 1 95 12 8 + time model var lon lat + 95 4 1 12 8 ``` + A new dimension will be added by running the function `InsertDim()` in order to obtain the same dimensions as the projections data. `InsertDim()` repeats the original data the required number of times (21 years of future simulations) in the adequated position (the temporal dimension in the summer_projection data is in the third position). ```r -climatology <- InsertDim(InsertDim(climatology, posdim = 2, lendim = 1), - posdim = 3, lendim = 95) +climatology <- InsertDim( + InsertDim( + climatology, posdim = 2, lendim = 1, name = 'var'), + posdim = 1, lendim = 95, name = 'time') ``` + The anomaly for each model is obtained by simply subtracting. @@ -214,7 +216,7 @@ anomaly <- summer_projection - climatology In order to obtain a spatial visualitzation, the temporal mean is computed. So, the time average anomalies for all models is saved in the `average` object. `AnoAgree()` function from **ClimProjDiags package** calculates the percentages of models which agrees with a positive or negative mean in each grid point. ```r -average <- Mean1Dim(anomaly, which(names(dim(anomaly)) == "time")) +average <- MeanDims(anomaly, dims = "time") agreement <- AnoAgree(average, membersdim = which(names(dim(average)) == "model")) ``` @@ -227,7 +229,7 @@ agreement_threshold <- 80 colorbar_lim <- ceiling(max(abs(max(average)), abs(min(average)))) brks <- seq(-colorbar_lim, colorbar_lim, length.out = 11) -PlotEquiMap(drop(Mean1Dim(average, which(names(dim(average)) == "model"))), +PlotEquiMap(drop(MeanDims(average, dims = "model")), lat = lat, lon = lon, units = "K", brks = brks, toptitle = paste(var, "- climatology:", start_climatology, "to", end_climatology, "and future simulation:", @@ -243,17 +245,17 @@ PlotEquiMap(drop(Mean1Dim(average, which(names(dim(average)) == "model"))), ### 5- Multi-model agreement temporal visualization -To visualize the time evolution of multi-model agreement, the spatial average is performed by a grid pixel size using the `WeightedMean` function from the **ClimProjDiags package**. Also, a smooth filter is applied with the `Smoothing()` function from the **s2dverifiction package**. In this example, a 5-year moving window filter is applied by defining the parameter `runmeanlen = 5`. +To visualize the time evolution of multi-model agreement, the spatial average is performed by a grid pixel size using the `WeightedMean` function from the **ClimProjDiags package**. Also, a smooth filter is applied with the `Smoothing()` function from the **s2dv package**. In this example, a 5-year moving window filter is applied by defining the parameter `runmeanlen = 5`. ```r temporal <- drop(WeightedMean(anomaly, lon = lon, lat = lat, mask = NULL)) -temporal <- Smoothing(temporal, runmeanlen = 5, numdimt = 2) +temporal <- Smoothing(temporal, time_dim = 'time', runmeanlen = 5) ``` Before visualizing, a data frame with the proper format is created. ```r -data_frame <- as.data.frame.table(t(temporal)) +data_frame <- as.data.frame.table(temporal) years <- rep(start_projection : end_projection, 4) data_frame$Year <- c(years) names(data_frame)[2] <- "Model" diff --git a/vignettes/diurnaltemp.Rmd b/vignettes/diurnaltemp.Rmd index 2abf5b5..e38963f 100644 --- a/vignettes/diurnaltemp.Rmd +++ b/vignettes/diurnaltemp.Rmd @@ -32,7 +32,7 @@ library(ClimProjDiags) All the other R packages involved can be installed directly from CRAN and loaded as follows: ```r -library(s2dverification) +library(s2dv) library(abind) library(parallel) ``` @@ -78,12 +78,16 @@ for (j in 1 : 8) { tmin_historical <- abind(tmin_historical, gridnew2, along = 3) } -tmax_historical <- InsertDim(InsertDim(tmax_historical, posdim = 1, lendim = 1), - posdim = 1, lendim = 1) -tmin_historical <- InsertDim(InsertDim(tmin_historical, posdim = 1, lendim = 1), - posdim = 1, lendim = 1) -names(dim(tmax_historical)) <- c("model", "var", "time", "lon", "lat") -names(dim(tmin_historical)) <- c("model", "var", "time", "lon", "lat") +names(dim(tmax_historical)) <- c('time', 'lon', 'lat') +names(dim(tmin_historical)) <- c('time', 'lon', 'lat') + +tmax_historical <- InsertDim(InsertDim(tmax_historical, posdim = 1, + lendim = 1, name = 'var'), + posdim = 1, lendim = 1, name = 'model') +tmin_historical <- InsertDim(InsertDim(tmin_historical, posdim = 1, + lendim = 1, name = 'var'), + posdim = 1, lendim = 1, name = 'model') + time <- seq(ISOdate(1971, 1, 1), ISOdate(2000, 12, 31), "day") metadata <- list(time = list(standard_name = 'time', long_name = 'time', calendar = 'proleptic_gregorian', @@ -125,12 +129,16 @@ for (j in 1 : 8) { tmin_projection <- abind(tmin_projection, gridnew2, along = 3) } -tmax_projection <- InsertDim(InsertDim(tmax_projection, posdim = 1, lendim = 1), - posdim = 1, lendim = 1) -tmin_projection <- InsertDim(InsertDim(tmin_projection, posdim = 1, lendim = 1), - posdim = 1, lendim = 1) -names(dim(tmax_projection)) <- c("model", "var", "time", "lon", "lat") -names(dim(tmin_projection)) <- c("model", "var", "time", "lon", "lat") +names(dim(tmax_projection)) <- c("time", "lon", "lat") +names(dim(tmin_projection)) <- c("time", "lon", "lat") + +tmax_projection <- InsertDim(InsertDim(tmax_projection, posdim = 1, + lendim = 1, name = 'var'), + posdim = 1, lendim = 1, name = 'model') +tmin_projection <- InsertDim(InsertDim(tmin_projection, posdim = 1, + lendim = 1, name = 'var'), + posdim = 1, lendim = 1, name = 'model') + time <- seq(ISOdate(2006, 1, 1), ISOdate(2100, 12, 31), "day") metadata <- list(time = list(standard_name = 'time', long_name = 'time', calendar = 'proleptic_gregorian', @@ -193,7 +201,7 @@ A four panel plot can be generated to visualize the seasonal indicator by runnin ```r dtr_rcp <- array(dim = c(length(lon), length(lat), 4)) for (j in 1:4){ - dtr_rcp[,,j] <- Mean1Dim(dtr_indicator$indicator[,j,,,,], 1) + dtr_rcp[,,j] <- MeanDims(dtr_indicator$indicator[,j,,,,], 'year') } breaks <- 0 : (max(dtr_rcp) + 1) @@ -224,10 +232,10 @@ dtr_indicator_reference <- DTRIndicator(tmax = tmax_historical, The comparison between the reference and the future projection diurnal temperature variation indicator can be computed With a simple subtraction. To visualize the result, `PlotLayout` function will be applied again. The resulting plot will be saved in the working directory. ```r -dtr_diff <- array(dim=c(length(lat), length(lon), 4)) +dtr_diff <- array(dim = c(length(lat), length(lon), 4)) for (i in 1:4){ - dtr_diff[,,i] <- Mean1Dim(dtr_indicator$indicator[,i,1,1,,], 1) - - Mean1Dim(dtr_indicator_reference$indicator[,i,1,1,,], 1) + dtr_diff[,,i] <- MeanDims(dtr_indicator$indicator[, i, 1, 1, , ], 'year') - + MeanDims(dtr_indicator_reference$indicator[, i, 1, 1, , ], 'year') } PlotLayout(PlotEquiMap, c(1, 2), lon = lon, lat = lat, var = dtr_diff, diff --git a/vignettes/extreme_indices.Rmd b/vignettes/extreme_indices.Rmd index 57eb34f..cbbbf12 100644 --- a/vignettes/extreme_indices.Rmd +++ b/vignettes/extreme_indices.Rmd @@ -39,7 +39,7 @@ All the other R packages involved can be installed directly from CRAN and loaded ```r -library(s2dverification) +library(s2dv) library(abind) library(multiApply) library(ggplot2) @@ -71,9 +71,11 @@ for (j in 1 : 8) { gridnew <- apply(gridlon, 2, function(x) {x - rnorm(10958, mean = j * 0.5, sd = 3)}) tmax_historical <- abind(tmax_historical, gridnew, along = 3) } -tmax_historical <- InsertDim(InsertDim(tmax_historical, posdim = 1, lendim = 1), - posdim = 1, lendim = 1) -names(dim(tmax_historical)) <- c("model", "var", "time", "lon", "lat") +names(dim(tmax_historical)) <- c("time", "lon", "lat") + +tmax_historical <- InsertDim(InsertDim(tmax_historical, posdim = 1, + lendim = 1, name = 'var'), + posdim = 1, lendim = 1, name = 'model') time <- seq(ISOdate(1971, 1, 1), ISOdate(2000, 12, 31), "day") metadata <- list(time = list(standard_name = 'time', long_name = 'time', calendar = 'proleptic_gregorian', @@ -102,8 +104,10 @@ for (j in 1 : 8) { rnorm(34698, mean = j * 0.5, sd = 3)}) tmax_projection <- abind(tmax_projection, gridnew, along = 3) } -tmax_projection <- InsertDim(InsertDim(tmax_projection, posdim = 1, lendim = 1), posdim = 1, lendim = 1) -names(dim(tmax_projection)) <- c("model", "var", "time", "lon", "lat") +names(dim(tmax_projection)) <- c("time", "lon", "lat") +tmax_projection <- InsertDim(InsertDim(tmax_projection, posdim = 1, + lendim = 1, name = 'var'), + posdim = 1, lendim = 1, name = 'model') time <- seq(ISOdate(2006, 1, 1), ISOdate(2100, 12, 31), "day") metadata <- list(time = list(standard_name = 'time', long_name = 'time', calendar = 'proleptic_gregorian', @@ -159,12 +163,12 @@ names(dim(anomaly_data))[1] <- "time" ``` -This data can be detrended by applying the `Trend` function from **s2dverification package**. In order to remove the trend from the `tmax_historical`, the correction is calculated by subtracting the `detrended_data` to the `anomaly_data`. +This data can be detrended by applying the `Trend` function from **s2dv package**. In order to remove the trend from the `tmax_historical`, the correction is calculated by subtracting the `detrended_data` to the `anomaly_data`. ```r -detrended_data <- Trend(anomaly_data, - posTR = which(names(dim(anomaly_data)) == "time")) +detrended_data <- Trend(anomaly_data, time_dim = "time") + diff <- anomaly_data - detrended_data$detrended diff <- aperm(diff, c(2,3,1,4,5)) detrended_data <- tmax_historical - diff @@ -247,8 +251,7 @@ masc <- rep(0, 8 * 12) masc[c(5 : 12, 18 : 24, 31 : 34, 43, 44, 47, 56 : 60, 67 : 72, 79, 82 : 84, 93 : 96)] <- 1 dim(masc) <- c(12, 8) -PlotEquiMap(Mean1Dim(HeatExtremeIndex, - which(names(dim(HeatExtremeIndex)) == "year")), +PlotEquiMap(MeanDims(HeatExtremeIndex, "year"), lon = lon, lat = lat, filled.continents = FALSE, toptitle = "Extreme Heat Index", dots = masc, fileout = "SpatialExtremeHeatIndex.png") @@ -263,7 +266,7 @@ The inland average of the Extreme Heat Index can be computed to plot its time ev ```r temporal <- WeightedMean(HeatExtremeIndex, lon = lon, lat = lat, mask = drop(masc)) -temporal_3ysmooth <- Smoothing(temporal, runmeanlen = 3, numdimt = 1) +temporal_3ysmooth <- Smoothing(temporal, runmeanlen = 3, time_dim = 'year') ``` @@ -335,8 +338,7 @@ Spatial representation of the Extreme Drought Index: ```r -PlotEquiMap(Mean1Dim(DroughtExtremeIndex, - which(names(dim(DroughtExtremeIndex)) == "year")), +PlotEquiMap(MeanDims(DroughtExtremeIndex, "year"), lon = lon, lat = lat, filled.continents = FALSE, toptitle = "Drought Index", brks = seq(-1, 1, 0.01), fileout = "SpatialDroughtIndex.png") @@ -352,7 +354,7 @@ Evolution of inland average of the Extreme Drought Index: ```r temporal <- WeightedMean(DroughtExtremeIndex, lon = lon, lat = lat, mask = drop(masc)) -temporal_5ysmooth <- Smoothing(temporal, runmeanlen = 5, numdimt = 1) +temporal_5ysmooth <- Smoothing(temporal, runmeanlen = 5, time_dim = "year") png("Temporal_Inland_ExtremeDroughtIndex.png", width = 8, height = 5, units= 'in', res = 100, type = "cairo") plot(2006: 2100, temporal, type = "l", lty = 5, lwd = 2, bty = 'n', @@ -402,9 +404,8 @@ Spatial representation of the Extreme Flooding Index: ```r -PlotEquiMap(Mean1Dim(FloodingExtremeIndex, - which(names(dim(FloodingExtremeIndex)) == "year")), lon = lon, - lat = lat, filled.continents = FALSE, +PlotEquiMap(MeanDims(FloodingExtremeIndex, "year"), + lon = lon, lat = lat, filled.continents = FALSE, toptitle = "Extreme Flooding Index", brks = seq(-1, 1, 0.1), fileout = "SpatialFloodingIndex.png") ``` @@ -419,7 +420,7 @@ PlotEquiMap(Mean1Dim(FloodingExtremeIndex, ```r temporal <- WeightedMean(FloodingExtremeIndex, lon = lon, lat = lat, mask = drop(masc)) -temporal_3ysmooth <- Smoothing(temporal, runmeanlen = 3, numdimt = 1) +temporal_3ysmooth <- Smoothing(temporal, runmeanlen = 3, time_dim = "year") png("Temporal_Inland_ExtremeFloodingIndex.png", width = 8, height = 5, units= 'in', res = 100, type = "cairo") plot(2006 : 2100, temporal, type = "l", lty = 5, lwd = 2, bty = 'n', @@ -460,7 +461,7 @@ A spatial visualization can be performed by executing: ```r -PlotEquiMap(Mean1Dim(aci, which(names(dim(aci)) == "year")), lon = lon, +PlotEquiMap(MeanDims(aci, "year"), lon = lon, lat = lat, filled.continents = FALSE, toptitle = "Indices Combination", fileout = "CombinedIndices.png") ``` diff --git a/vignettes/heatcoldwaves.Rmd b/vignettes/heatcoldwaves.Rmd index c50b4c4..01c1ac5 100644 --- a/vignettes/heatcoldwaves.Rmd +++ b/vignettes/heatcoldwaves.Rmd @@ -32,7 +32,7 @@ All the other R packages involved can be installed directly from CRAN and loaded ```r library(abind) -library(s2dverification) +library(s2dv) library(parallel) ``` @@ -68,9 +68,10 @@ for (j in 1 : 8) { gridnew <- apply(gridlon, 2, function(x) {x - rnorm(10958, mean = j * 0.5, sd = 3)}) tmax_historical <- abind(tmax_historical, gridnew, along = 3) } -tmax_historical <- InsertDim(InsertDim(tmax_historical, posdim = 1, lendim = 1), - posdim = 1, lendim = 1) -names(dim(tmax_historical)) <- c("model", "var", "time", "lon", "lat") +names(dim(tmax_historical)) <- c("time", "lon", "lat") +tmax_historical <- InsertDim(InsertDim(tmax_historical, posdim = 1, + lendim = 1, name = "var"), + posdim = 1, lendim = 1, name = "model") time <- seq(ISOdate(1971, 1, 1), ISOdate(2000, 12, 31), "day") metadata <- list(time = list(standard_name = 'time', long_name = 'time', calendar = 'proleptic_gregorian', @@ -97,9 +98,10 @@ for (j in 1 : 8) { sd = 3)}) tmax_projection <- abind(tmax_projection, gridnew, along = 3) } -tmax_projection <- InsertDim(InsertDim(tmax_projection, posdim = 1, lendim = 1), - posdim = 1, lendim = 1) -names(dim(tmax_projection)) <- c("model", "var", "time", "lon", "lat") +names(dim(tmax_projection)) <- c("time", "lon", "lat") +tmax_projection <- InsertDim(InsertDim(tmax_projection, posdim = 1, + lendim = 1, name = "var"), + posdim = 1, lendim = 1, name = "model") time <- seq(ISOdate(2006, 1, 1), ISOdate(2100, 12, 31), "day") metadata <- list(time = list(standard_name = 'time', long_name = 'time', calendar = 'proleptic_gregorian', -- GitLab