feat: refactor SoilGrids functions to handle multiple properties #3455
Conversation
…improve efficiency
|
I made changes to refactor SoilGrids functions for improved efficiency and scalability; please take a look,i think this works for you. |
|
@DongchenZ @Qianyuxuan Can you review this PR and I think it solves the mentioned issue.Just i want you to go through it once and let me know if you feel any changes required. |
There was a problem hiding this comment.
This PR improves the ability to extract soil properties other than SOC. However, I worry that those functions (which heavily rely on the terra package) will be limited to broader spatial extents (e.g., North America, Asia, Africa, etc.) where the computation demand and the memory usage will be unrealistic high. Given that, I think it will be helpful to have a prior check on the volume of data requested (e.g., if the size of the data is too large, then never try it and just return an error message).
For larger spatial extents, it would be more efficient to download the tiled SoilGrids data directly (I usually use the wget command) and merge the tiled images using the gdalwarp command.
| #' @param local_raster Optional. The file path of a local raster to query. If NULL, queries the remote VRT. | ||
| #' @return A data frame containing the queried data | ||
| #' @export | ||
| query_soilgrids <- function(points, soil_property, depth, quantile, local_raster = NULL) { |
There was a problem hiding this comment.
See my previous comments for projecting the points from lat/lon degrees to the SoilGrids coordinate system.
| dplyr::summarize( | ||
| mean = mean(value, na.rm = TRUE), | ||
| sd = sd(value, na.rm = TRUE), | ||
| quantile_05 = quantile(value, probs = quantiles[1], na.rm = TRUE), |
There was a problem hiding this comment.
I would not have those hardcoded names for different quantiles. Instead, you might consider calculating them first and naming them afterward based on the quantile we choose.
There was a problem hiding this comment.
yes, I'll adjust the code to calculate and name the quantiles dynamically
| #' @param quantiles A vector of quantiles to estimate uncertainties (default: c(0.05, 0.5, 0.95)) | ||
| #' @return A data frame containing the uncertainties | ||
| #' @export | ||
| soilgrids_uncertainty <- function(data, quantiles = c(0.05, 0.5, 0.95)) { |
There was a problem hiding this comment.
I wonder if the calculation of uncertainty here is correct or not. Our current calculations will first simulate the uncertainty of soil properties across the vertical profile using the mean and different quantiles from SoilGrids directly and then summarize the uncertainty once we have a good fit on that (See
There was a problem hiding this comment.
I think @Qianyuxuan might have better opinions on that.
There was a problem hiding this comment.
I am confused about how you calculate sd here: "sd = sd(value, na.rm = TRUE)" with known quantiles and means. Usually, we estimate means and sds based on the best fits. For organic carbon density (ocd), it was found to follow gamma distribution best based on reported means and quantiles. But for soil texture data (sand, clay, and silt fractions), they follow Dirichlet distribution and estimate of corresponding parameters can be found in my pending PR here: https://github.com/PecanProject/pecan/pull/3406/files#diff-8b26bae8174687f7e2bf631c5188b979a5b8e4bdacc5c604da84e09a2abc0058
There was a problem hiding this comment.
If I'm reading the current code correctly, it's grouping by site, and thus calculating a mean and sd over depth, which is wrong on many levels (pun intended). Sampling from rnorm is also wrong. This code is in no way a "refactor" of the existing code
…ds_carbonstocks function
That's an important point! Implementing a data volume check could help manage resource usage effectively. Your approach of using tiled downloads and merging with |
| upper_bound <- segment[2] | ||
|
|
||
| total_soc <- soc_data %>% | ||
| dplyr::filter(depths >= lower_bound & depths < upper_bound) %>% |
There was a problem hiding this comment.
Is “soc_data” soil organic carbon content (SOC, in g/kg) or organic carbon density (OCD, in kg/m³)? I guess it is OCD if you multiply it by thickness (see https://www.isric.org/explore/soilgrids/faq-soilgrids). For total_SOC, it should be something like "total_SOC=sum(OCDthickness(1-coarse_fraction))" but can you explain your calculation here?
There was a problem hiding this comment.
The soc_data refers to organic carbon density (OCD) in kg/m³. The calculation for total_SOC is performed as follows:
[
\text{total_SOC} = \sum(\text{OCD} \times \text{thickness} \times (1 - \text{coarse_fraction}))
]
This accounts for the volume of soil represented by the thickness and adjusts for the coarse fraction to provide an accurate estimate of total SOC.
There was a problem hiding this comment.
- your math seems right, but that math isn't what's implemented in the code
- I'm not sure the "if" will work as intended for vector data
- The "if" assumes that if the coarse fraction is NULL then it's 0. Setting aside that this might not be a valid assumption, it would be far simpler to do that substitution earlier in the pipe
|
@Qianyuxuan @DongchenZ any insights on this changes? |
AritraDey-Dev
changed the title
4 participants
fixes #3416
Refactored the SoilGrids functions in
modules/data.land/R/soilgrids_utils.Rto enhance efficiency and scalability. The changes include:download_soilgrids_rasterto handle multiple soil properties.query_soilgridsto query data for multiple soil properties.soilgrids_carbonstocksto include an optional coarse fraction in carbon stock calculations.