1- # #' A function to estimate individual alphas for Dirichlet distribution to approximate the observed quantiles with means as known moments
2- # #' for SoilGrids soil texture data.
3- # #' Dirichlet distribution is assumed as soil texture data follow categorical distribution and the probability of each category is in the range 0 to 1,
4- # #' and all must sum to 1.
5- # #'
6- # #' @param means A vector of means of sand, clay, and silt proportion for one soil layer at one site from SoilGrids data
7- # #' @param quantiles A list of 5th, 50th, and 95th percentiles for sand, clay, and silt for one soil layer at one site from SoilGrids data
8- # #'
9- # #' @example
10- # #' \dontrun{
11- # #' means <- c(0.566,0.193,0.241) # means of sand, clay, and silt at one site and one depth
12- # #' quantiles <-list(
13- # #' q5 = c(0.127,0.034,0.052), # 5th percentile for each category: sand, clay, and silt at one site and one depth
14- # #' q50 = c(0.615,0.15,0.191), # 50th percentile (median) for each category: sand, clay, and silt at one site and one depth
15- # #' q95 = c(0.799,0.66,0.616)) # 95th percentile for each category: sand, clay, and silt at one site and one depth
16- # #' alpha_est <- estimate_dirichlet_parameters(means, quantiles)
17- # #' }
18- # #' @return The individual alphas that work best to fit the observed quantiles
19- # #' @export
20- # #' @author Qianyu Li
1+ # ' A function to estimate individual alphas for Dirichlet distribution to
2+ # ' approximate the observed quantiles with means as known moments for SoilGrids
3+ # ' soil texture data.
4+ # ' Dirichlet distribution is assumed as soil texture data follow categorical
5+ # ' distribution and the probability of each category is in the range 0 to 1,
6+ # ' and all must sum to 1.
7+ # '
8+ # ' @param means A vector of means of sand, clay, and silt proportion for one
9+ # ' soil layer at one site from SoilGrids data
10+ # ' @param quantiles A list of 5th, 50th, and 95th percentiles for sand, clay,
11+ # ' and silt for one soil layer at one site from SoilGrids data
12+ # '
13+ # ' @examples
14+ # ' \dontrun{
15+ # ' # Means and percentiles for each category: sand, clay, and silt at one site and one depth
16+ # ' means <- c(0.566,0.193,0.241)
17+ # ' quantiles <-list(
18+ # ' q5 = c(0.127,0.034,0.052), # 5th percentile
19+ # ' q50 = c(0.615,0.15,0.191), # 50th percentile (median)
20+ # ' q95 = c(0.799,0.66,0.616)) # 95th percentile
21+ # ' alpha_est <- estimate_dirichlet_parameters(means, quantiles)
22+ # ' }
23+ # ' @return The individual alphas that work best to fit the observed quantiles
24+ # ' @author Qianyu Li
2125estimate_dirichlet_parameters <- function (means , quantiles ) {
2226
2327 # A function to optimize alpha0, which is the sum of individual alphas.
@@ -31,13 +35,19 @@ estimate_dirichlet_parameters <- function(means, quantiles) {
3135 # Generate samples based on estimated alpha
3236 samples <- MCMCpack :: rdirichlet(10000 , alpha ) # Generate samples
3337 # Compute differences with observed quantiles
34- estimated_quantiles <- apply(samples , 2 , quantile , probs = c(0.05 , 0.5 , 0.95 ),na.rm = TRUE )
38+ estimated_quantiles <- apply(
39+ x = samples ,
40+ margin = 2 ,
41+ FUN = stats :: quantile ,
42+ probs = c(0.05 , 0.5 , 0.95 ),
43+ na.rm = TRUE
44+ )
3545 quantile_diff <- sum((estimated_quantiles - do.call(rbind , quantiles ))^ 2 )
3646 return (quantile_diff )
3747 }
3848
3949 # Optimize alpha0
40- result <- optim(
50+ result <- stats :: optim(
4151 par = 1 , # Initial guess for alpha0
4252 fn = objective_function ,
4353 method = " L-BFGS-B" ,
@@ -56,31 +66,34 @@ estimate_dirichlet_parameters <- function(means, quantiles) {
5666
5767
5868
59- # #' A function to estimate the soil parameters based on SoilGrids soil texture data and write the parameter paths into settings
60- # #'
61- # #' @param settings A multi-site settings
62- # #' @param sand A data frame containing sand fraction in percentage from SoilGrids250m v2.0 with columns "Depth", "Quantile", "Siteid", and "Value"
63- # #' @param clay A data frame containing clay fraction in percentage from SoilGrids250m v2.0 with columns "Depth", "Quantile", "Siteid", and "Value"
64- # #' @param silt A data frame containing silt fraction in percentage from SoilGrids250m v2.0 with columns "Depth", "Quantile", "Siteid", and "Value"
65- # #' @param outdir Provide the path to store the parameter files
66- # #' @param write_into_settings Whether to write the path of parameter file into the setting. The default is TRUE
67- # #'
68- # #' @example
69- # #' \dontrun{
70- # #'
71- # #' outdir <- "/projectnb/dietzelab/Cherry/SoilGrids_texture/39NEON"
72- # #' sand <- readRDS("/projectnb/dietzelab/Cherry/SoilGrids_texture/sand_percent.rds") #sand fraction in percentage
73- # #' clay <- readRDS("/projectnb/dietzelab/Cherry/SoilGrids_texture/clay_percent.rds") #clay fraction in percentage
74- # #' silt <- readRDS("/projectnb/dietzelab/Cherry/SoilGrids_texture/silt_percent.rds") #silt fraction in percentage
75- # #' settings <-read.settings("/projectnb/dietzelab/Cherry/xml/pecan_monthly_SDA_soilwater.xml")
76- # #' soil_params_ensemble_soilgrids(settings,sand,clay,silt,outdir)
77- # #' }
78- # #'
79- # #' @return Ensemble soil parameter files defined in outdir and file paths in xml file
80- # #' @export
81- # #' @author Qianyu Li
82- # #' @importFrom magrittr %>%
83- # #'
69+ # ' A function to estimate the soil parameters based on SoilGrids soil texture
70+ # ' data and write the parameter paths into settings
71+ # '
72+ # ' @param settings A multi-site settings
73+ # ' @param sand,clay,silt Data frames containing fraction in percentage from SoilGrids250m
74+ # ' v2.0, each with columns "Depth", "Quantile", "Siteid", and "Value"
75+ # ' @param outdir Provide the path to store the parameter files
76+ # ' @param write_into_settings Whether to write the path of parameter file into
77+ # ' the setting. The default is TRUE
78+ # '
79+ # ' @examples
80+ # ' \dontrun{
81+ # '
82+ # ' outdir <- "/projectnb/dietzelab/Cherry/SoilGrids_texture/39NEON"
83+ # ' # each file contains percent salt, silt, or clay
84+ # ' sand <- readRDS("/path/to/SoilGrids_texture/sand_percent.rds")
85+ # ' clay <- readRDS("/path/to/SoilGrids_texture/clay_percent.rds")
86+ # ' silt <- readRDS("/path/to/SoilGrids_texture/silt_percent.rds")
87+ # ' settings <-read.settings("/path/to/pecan_monthly_SDA_soilwater.xml")
88+ # ' soil_params_ensemble_soilgrids(settings,sand,clay,silt,outdir)
89+ # ' }
90+ # '
91+ # ' @return Ensemble soil parameter files defined in outdir and file paths in xml file
92+ # ' @export
93+ # ' @author Qianyu Li
94+ # ' @importFrom magrittr %>%
95+ # ' @importFrom rlang .data
96+ # ' @importFrom foreach %dopar%
8497
8598soil_params_ensemble_soilgrids <- function (settings ,sand ,clay ,silt ,outdir ,write_into_settings = TRUE ){
8699
@@ -154,21 +167,28 @@ soil_params_ensemble_soilgrids <- function(settings,sand,clay,silt,outdir,write_
154167 # Estimate Dirichlet parameters for each depth at each site
155168 for (depths in sort(unique(texture_all $ soil_depth ))) {
156169 quantiles <- list (
157- q5 = dplyr :: filter(dat [[i ]], quantile == " 0.05" , soil_depth == depths ) %> % dplyr :: select(
158- fraction_of_sand_in_soil ,
159- fraction_of_clay_in_soil ,
160- fraction_of_silt_in_soil ), # 5th percentile for each category
161- q50 = dplyr :: filter(dat [[i ]], quantile == " 0.5" , soil_depth == depths ) %> % dplyr :: select(
162- fraction_of_sand_in_soil ,
163- fraction_of_clay_in_soil ,
164- fraction_of_silt_in_soil ), # 50th percentile (median) for each category
165- q95 = dplyr :: filter(dat [[i ]], quantile == " 0.95" , soil_depth == depths ) %> % dplyr :: select(
166- fraction_of_sand_in_soil ,
167- fraction_of_clay_in_soil ,
168- fraction_of_silt_in_soil )) # 95th percentile for each category
170+ q5 = dplyr :: filter(dat [[i ]], .data $ quantile == " 0.05" , soil_depth == depths ) %> %
171+ dplyr :: select(
172+ " fraction_of_sand_in_soil" ,
173+ " fraction_of_clay_in_soil" ,
174+ " fraction_of_silt_in_soil" ), # 5th percentile for each category
175+ q50 = dplyr :: filter(dat [[i ]], .data $ quantile == " 0.5" , soil_depth == depths ) %> %
176+ dplyr :: select(
177+ " fraction_of_sand_in_soil" ,
178+ " fraction_of_clay_in_soil" ,
179+ " fraction_of_silt_in_soil" ), # 50th percentile (median) for each category
180+ q95 = dplyr :: filter(dat [[i ]], .data $ quantile == " 0.95" , soil_depth == depths ) %> %
181+ dplyr :: select(
182+ " fraction_of_sand_in_soil" ,
183+ " fraction_of_clay_in_soil" ,
184+ " fraction_of_silt_in_soil" )) # 95th percentile for each category
169185
170186 # Extract the means
171- means <- dplyr :: filter(dat [[i ]], quantile == " Mean" , soil_depth == depths ) %> % dplyr :: select(fraction_of_sand_in_soil ,fraction_of_clay_in_soil ,fraction_of_silt_in_soil )
187+ means <- dplyr :: filter(dat [[i ]], .data $ quantile == " Mean" , soil_depth == depths ) %> %
188+ dplyr :: select(
189+ " fraction_of_sand_in_soil" ,
190+ " fraction_of_clay_in_soil" ,
191+ " fraction_of_silt_in_soil" )
172192 soil_rescaled <- rescale_sum_to_one(means $ fraction_of_sand_in_soil ,means $ fraction_of_clay_in_soil ,means $ fraction_of_silt_in_soil )
173193
174194 # Replace the original means with the rescaled ones
@@ -182,7 +202,7 @@ soil_params_ensemble_soilgrids <- function(settings,sand,clay,silt,outdir,write_
182202 # Generate the ensemble soil texture data based on the ensemble size (ens_n) defined in the settings
183203 samples <- MCMCpack :: rdirichlet(ens_n , alpha_est )
184204 colnames(samples ) <- c(" fraction_of_sand_in_soil" ," fraction_of_clay_in_soil" ," fraction_of_silt_in_soil" )
185- samples <- list (samples ) %> % setNames(depths )
205+ samples <- list (samples ) %> % stats :: setNames(depths )
186206 samples_ens <- append(samples_ens , samples )
187207 }
188208 # Generate soil parameter file for each one in ensemble soil texture data
@@ -211,7 +231,7 @@ soil_params_ensemble_soilgrids <- function(settings,sand,clay,silt,outdir,write_
211231 settings [[ind ]]$ run $ inputs $ soil_physics $ ensemble <- ens_n
212232 settings [[ind ]]$ run $ inputs $ soil_physics $ path <- create_mult_list(rep(" path" , ens_n ), PATH [[i ]])
213233 }
214- write.settings(settings ,outputdir = settings $ outdir ,outputfile = " pecan.xml" )
234+ PEcAn.settings :: write.settings(settings ,outputdir = settings $ outdir ,outputfile = " pecan.xml" )
215235 }
216236}
217237
@@ -224,7 +244,7 @@ reformat_soil_list <- function(samples_all_depth) {
224244 " fraction_of_silt_in_soil" )
225245
226246 # Initialize a new list to store reformatted data
227- reformatted <- setNames(vector(" list" , length(fractions )), fractions )
247+ reformatted <- stats :: setNames(vector(" list" , length(fractions )), fractions )
228248
229249 # Extract data for each fraction
230250 for (fraction in fractions ) {
0 commit comments