diff --git a/inst/doc/faq.md b/inst/doc/faq.md index ea84377aa7425394a66bef20794a39072c9451c5..689e13a2f468cc1ccf27ffb146af46047794a2c4 100644 --- a/inst/doc/faq.md +++ b/inst/doc/faq.md @@ -37,7 +37,7 @@ This document intends to be the first reference for any doubts that you may have 4. [My jobs work well in workstation and fatnodes but not on Power9 (or vice versa)](#4-my-jobs-work-well-in-workstation-and-fatnodes-but-not-on-power9-or-vice-versa) 5. [Errors related to wrong file formatting](#5-errors-related-to-wrong-file-formatting) 6. [Errors using a new cluster (setting Nord3)](#6-errors-using-a-new-cluster-setting-nord3) - 7. [Start() fails retrieving data](#7-start-fails-retrieving-data) + 7. [Start() fails retrieving data](#7-start-fails-retrieving-data) 8. [Error 'caught segfault' when job submitted to HPCs](#8-error-caught-segfault-when-job-submitted-to-hpcs) @@ -1123,3 +1123,13 @@ the process you are running may fail. Remove the files and re-run your code. If the error persists, check your code with a smaller data sample to discard a problem with your code since this error message indicates that you are requesting more memory than the available. + +### 9. White lines on the figure of interpolated irregular data + +To process irregular grid data, we can load the data by Start() and interpolate it to a regular grid by other tools (e.g., using s2dv::CDORemap). +In some instances, when we plot the interpolated data, we see some white lines on the map (see figure below). +To solve this problem, we can try to exclude the first and last indices of the latitude and longitude in Start() call. +Check [use case ex1_15](inst/doc/usecase/ex1_15_irregular_grid_CDORemap.R) for the example script (Case 2). The solution of each case may differ, so if you find this solution does not work for your case, please open an issue. + + + diff --git a/inst/doc/figures/faq_2_9_white_stripes.png b/inst/doc/figures/faq_2_9_white_stripes.png new file mode 100644 index 0000000000000000000000000000000000000000..ceccedbf1cefc0f8729425bae8c0b70a55cc4748 Binary files /dev/null and b/inst/doc/figures/faq_2_9_white_stripes.png differ diff --git a/inst/doc/usecase.md b/inst/doc/usecase.md index 9d1713901b754e36661493ac165d9bd7dc213466..0ffacc7f618336a7253c121a56a0c0f36dfdc303 100644 --- a/inst/doc/usecase.md +++ b/inst/doc/usecase.md @@ -65,6 +65,9 @@ parameter and Start() can recognize this dependecy relationship. the dependency between one inner dimension and one file dimension (i.e., the usage of *_across), while this use case is for two file dimensions (i.e., the usage of *_depends). + 15. [Load irregular grid data](inst/doc/usecase/ex1_15_irregular_grid_CDORemap.R) + This script shows how to use Start() to load irregular grid data , then regrid it by s2dv::CDORemap. + 2. **Execute computation (use `Compute()`)** 1. [Function working on time dimension](inst/doc/usecase/ex2_1_timedim.R) 2. [Function using attributes of the data](inst/doc/usecase/ex2_2_attr.R) @@ -74,23 +77,28 @@ this use case is for two file dimensions (i.e., the usage of *_depends). 4. [Use two functions in workflow](inst/doc/usecase/ex2_4_two_func.R) 5. 6. [Use external parameters in atomic function](inst/doc/usecase/ex2_6_ext_param_func.R) + 7. [Calculate the ensemble-adjusted Continuous Ranked Probability Score (CRPS)](inst/doc/usecase/ex2_7_seasonal_forecast_crps.R) Use `SpecsVerification::EnsCrps` to calculate the ensemble-adjusted Continuous Ranked Probability Score (CRPS) for ECWMF experimental data, and do ensemble mean. Use `s2dv::PlotEquiMap` to plot the CRPS map. + 8. [Use CSTools Calibration function](inst/doc/usecase/ex2_8_calibration.R) Use `CSTools:::.cal`, the interior function of `CSTools::CST_Calibration`, to do the bias adjustment for ECMWF experimental monthly mean data. + 9. [Use a mask to apply different methods to different gridpoints](inst/doc/usecase/ex2_9_mask.R) If you need to apply your analysis in a few gridpoints, you may want to consider use case 1.6, but if you need to load a lot of grid points, maybe this a better solution. + 10. [Apply an existing mask on data](inst/doc/usecase/ex2_10_existing_mask.R) This use case shows you how to apply the existing mask file on your data. If you need to create the mask file on your own, go to ex2_9_mask.R. + 11. [Two datasets with different length of target dimensions](inst/doc/usecase/ex2_11_two_dat_inconsistent_target_dim.R) This use case uses experimental and the corresponding observational data to calculate the temporal mean and spatial weighted mean. Notice that the spatial resolutions of the two datasets are different, but it still works because lat and lon are target dimensions. - 12. [Transform and chunk over spatial dimensions](inst/doc/usecase/ex2_12_transform_and_chunk.R) - This use case provides an example of transforming and chunking -latitude and longitude dimensions. If all other dimensions are used as target dimensions in the operation, -it would be good to have the option of chunking the spatial dimensions. - 13. [Interpolate irregular grid in the workflow](inst/doc/usecase/ex2_13_irregular_regrid.R) + + 12. [Transform and chunk spatial dimensions](inst/doc/usecase/ex2_12_transform_and_chunk.R) + This use case provides an example of transforming and chunking latitude and longitude dimensions. + + 13. [Load irregular grid data and interpolate it in the workflow](inst/doc/usecase/ex2_13_irregular_regrid.R) This script shows how to load irregular grid data by Start(), then regrid it -by s2dv::CDORemap in the workflow. It is a solution before Start() can deal with irregular regridding directly. +by s2dv::CDORemap in the workflow. It is a solution before Start() can deal with irregular regridding directly. diff --git a/inst/doc/usecase/ex1_15_irregular_grid_CDORemap.R b/inst/doc/usecase/ex1_15_irregular_grid_CDORemap.R new file mode 100644 index 0000000000000000000000000000000000000000..f8634430e13730ea135182d6b3da0206f96ee918 --- /dev/null +++ b/inst/doc/usecase/ex1_15_irregular_grid_CDORemap.R @@ -0,0 +1,83 @@ +#---------------------------------------------------------------------------- +# Author: An-Chi Ho +# Date: 18th Jun 2021 +# +# This script shows how to use Start() to load irregular grid data , then regrid it +# by s2dv::CDORemap. +#---------------------------------------------------------------------------- + +library(startR) +library(s2dv) + +# Case 1: +path1 <- '/esarchive/exp/ecearth/a1ua/cmorfiles/DCPP/EC-Earth-Consortium/EC-Earth3/dcppA-hindcast/r1i1p1f1/Omon/$var$/gn/v20190713/$var$_*_s$sdate$-$member$_gn_$aux$.nc' + +data1 <- Start(dataset = path1, + var = 'tos', + sdate = paste0(1960), + aux = indices(1), + aux_depends = 'sdate', + j = 'all', + i = 'all', + time = indices(1), + member = 'r1i1p1f1', + return_vars = list(j = NULL, i = NULL, + latitude = NULL, longitude = NULL), + retrieve = T) +lons1 <- attributes(data1)$Variables$common$longitude +lats1 <- attributes(data1)$Variables$common$latitude +dim(lons1) +# i j +#362 292 +dim(lats1) +# i j +#362 292 +dim(data1) +#dataset var sdate aux j i time member +# 1 1 1 1 292 362 1 1 + +res1 <- CDORemap(drop(data1), lons1, lats1, grid = 'r100x50', method = 'bil', crop = FALSE) +dim(res1$data_array) +#lon lat +#100 50 +PlotEquiMap(t(res1$data_array), lon = res1$lons, lat = res1$lats, fileout = '~/irre_grid_res1.png') + +#--------------------------------------------------------- + +# Case 2: +path2 <- paste0('/esarchive/exp/CMIP6/dcppA-hindcast/cmcc-cm2-sr5/cmip6-dcppA-hindcast_i1p1/', + 'DCPP/CMCC/CMCC-CM2-SR5/dcppA-hindcast/$member$/Omon/$var$/gn/v20210312/', + '$var$_*_s$sdate$-$member$_gn_$aux$.nc') + +data2 <- Start(dataset = path2, + var = 'tos', + sdate = '1960', + aux = 'all', + aux_depends = 'sdate', + j = indices(2:361), # remove two indices to avoid white strips + i = indices(2:291), # remove two indices to avoid white strips + time = indices(1), + member = 'r1i1p1f1', + return_vars = list(j = NULL, i = NULL, + latitude = NULL, longitude = NULL), + retrieve = T) + +lons2 <- attributes(data2)$Variables$common$longitude +lats2 <- attributes(data2)$Variables$common$latitude +dim(lons2) +# j i +#360 290 +dim(lats2) +# j i +#360 290 +dim(data2) +#dataset var sdate aux j i time member +# 1 1 1 1 360 290 1 1 + +data2 <- ClimProjDiags::Subset(data2, c(1,2,3,7,8), list(1,1,1,1,1), drop = 'selected') +res2 <- CDORemap(drop(data2), lons2, lats2, grid = 'r100x50', method = 'bil', crop = FALSE) +dim(res2$data_array) +#lon lat +#100 50 +PlotEquiMap(res2$data_array, lon = res2$lons, lat = res2$lats, fileout = '~/irre_grid_res2.png') + diff --git a/inst/doc/usecase/ex2_13_irregular_regrid.R b/inst/doc/usecase/ex2_13_irregular_regrid.R index df5e21f03aaeb3b47073d931f838e299334a1e69..e09a691333ad40754315d98e548fda55ba7aa9f0 100644 --- a/inst/doc/usecase/ex2_13_irregular_regrid.R +++ b/inst/doc/usecase/ex2_13_irregular_regrid.R @@ -8,6 +8,7 @@ #---------------------------------------------------------------------------- library(startR) +library(s2dv) path <- paste0('/esarchive/exp/CMIP6/dcppA-hindcast/cmcc-cm2-sr5/cmip6-dcppA-hindcast_i1p1/', 'DCPP/CMCC/CMCC-CM2-SR5/dcppA-hindcast/$member$/Omon/$var$/gn/v20210312/', @@ -26,6 +27,16 @@ data <- Start(dataset = path, latitude = NULL, longitude = NULL), retrieve = F) +attr(data, 'Dimensions') +#dataset var sdate aux j i time member +# 1 1 2 1 360 290 12 1 +dim(attr(data, 'Variables')$common$longitude) +# j i +#360 290 +dim(attr(data, 'Variables')$common$latitude) +# j i +#360 290 + func_regrid <- function(data) { lons <- attr(data, 'Variables')$common$longitude lats <- attr(data, 'Variables')$common$latitude @@ -36,6 +47,8 @@ func_regrid <- function(data) { return(list(data = data[[1]], lats = lats_reg, lons = lons_reg)) } +#NOTE: The data transposes if target_dims are only 'j' and 'i'. +# If only 'j' and 'i', output_dims will be 'lat', 'lon'. step <- Step(fun = func_regrid, target_dims = list(data = c('j', 'i')), output_dims = list(data = c('lon', 'lat'), @@ -58,11 +71,12 @@ dim(res$lats) # lat dataset var sdate aux time member # 180 1 1 2 1 12 1 -library(s2dv) -PlotEquiMap(drop(res$data)[ , , 1, 1], lon = drop(res$lons)[, 1, 1], +PlotEquiMap(drop(res$data)[ , , 1, 1], + lon = drop(res$lons)[, 1, 1], lat = drop(res$lats)[, 1, 1]) # Plot Layout for sdate = 1 all the time steps var <- Reorder(drop(res$data)[, , 1, ], c(3, 1, 2)) PlotLayout(PlotEquiMap, c('lon', 'lat'), var = var, lon = drop(res$lons)[, 1, 1], lat = drop(res$lats)[, 1, 1]) + diff --git a/tests/testthat/test-Compute-irregular_regrid.R b/tests/testthat/test-Compute-irregular_regrid.R new file mode 100644 index 0000000000000000000000000000000000000000..928fa493be57d7b9f975c57bbf229f986dab6d94 --- /dev/null +++ b/tests/testthat/test-Compute-irregular_regrid.R @@ -0,0 +1,71 @@ +library(s2dv) + +context("Irregular regriding in the workflow") + +test_that("1. ex2_13", { + +path <- paste0('/esarchive/exp/CMIP6/dcppA-hindcast/cmcc-cm2-sr5/cmip6-dcppA-hindcast_i1p1/', + 'DCPP/CMCC/CMCC-CM2-SR5/dcppA-hindcast/$member$/Omon/$var$/gn/v20210312/', + '$var$_*_s$sdate$-$member$_gn_$aux$.nc') +suppressWarnings( +data <- Start(dataset = path, + var = 'tos', + sdate = '1960', + aux = 'all', + aux_depends = 'sdate', + j = indices(2:361), # remove two indices to avoid white strips + i = indices(2:291), # remove two indices to avoid white strips + time = indices(1), + member = 'r1i1p1f1', + return_vars = list(j = NULL, i = NULL, + latitude = NULL, longitude = NULL), + retrieve = F) +) +func_regrid <- function(data) { + lons <- attr(data, 'Variables')$common$longitude + lats <- attr(data, 'Variables')$common$latitude + data <- s2dv::CDORemap(data, lons, lats, grid = 'r360x180', + method = 'bil', crop = FALSE) + lons_reg <- data[['lons']] + lats_reg <- data[['lats']] + return(list(data = data[[1]], lats = lats_reg, lons = lons_reg)) +} + +#NOTE: The data transposes if target_dims are only 'j' and 'i'. +# If only 'j' and 'i', output_dims will be 'lat', 'lon'. +step <- Step(fun = func_regrid, + target_dims = list(data = c('j', 'i')), + output_dims = list(data = c('lon', 'lat'), + lats = 'lat', lons = 'lon'), + use_attributes = list(data = "Variables")) +suppressWarnings( +wf <- AddStep(data, step) +) +suppressWarnings( +res <- Compute(workflow = wf$data, chunks = list(sdate = 1)) +) + +expect_equal( +dim(res$data), +c(lon = 360, lat = 180, dataset = 1, var = 1, sdate = 1, aux = 1, time = 1, member = 1) +) +expect_equal( +dim(res$lons), +c(lon = 360, dataset = 1, var = 1, sdate = 1, aux = 1, time = 1, member = 1) +) +expect_equal( +attr(data, 'Dimensions'), +c(dataset = 1, var = 1, sdate = 1, aux = 1, j = 360, i = 290, time = 1, member = 1) +) +expect_equal( +mean(res$data, na.rm = T), +13.20951, +tolerance = 0.0001 +) +expect_equal( +drop(res$data)[120,105:110], +c(28.32521, 28.07044, 27.59033, 27.02514, 26.55184, 26.67090), +tolerance = 0.0001 +) + +})