I know the issue of resampling has been raised before but I’ve noticed something which concerned me a bit in my latest use of the SDR model.
I am aware that input datasets are resampled to the resolution of the DEM, which in this case is 30 m. For the R factor, I have been using the global R Factor map from ESDAC which is recommended in the user guide and forums. As many of you will know gives rainfall erosivity in MJ mm ha-1 h-1 yr-1 at 1 km resolution.
My worries emerged when examining the “aligned_erosivity” raster intermediary output. I see the cell size for this raster was changed correctly to match the DEM, however the values on each pixel are the same as the original raster (aside from some interpolation). I could be getting this wrong, but in my mind should the per pixel value not go down when this dataset is resampled to 30 m resolution? For example, if the initial value is 1000 MJ mm/ha/year, should it not be divided by 11.1 (i.e. 100^2 / 30^2) in order to convert the per hectare value to a per pixel value (at 30 m resolution) of around 90 MJ mm/ha/year?
I am not sure whether my thinking is correct here but my worries emerged when I was comparing the USLE raster output from the model with other USLE studies. The InVEST output seemed like it was an order of magnitude or more too high. For example over Lesotho, the average USLE value from the model is 16.9 t/pixel/year (30 m), which translates to around 188 t/ha/year! This is confirmed by doing Zonal stats of my USLE layer for the whole country, which gives a figure of around 600 Million t of soil loss across Lesotho (which has an area of a little over 3 million ha). Lesotho is a mountainous country but this still seems really high. It’s almost as if the per pixel value is in fact the per ha value, and the total erosion is thus 11 times too high. My settings for the C and P factors seemed reasonable (I have added a file with a land cover legend and these values for reference). My K factor raster shouldn’t be the issue either, as the maximum value of this is around 0.05.
Given the error seems so big, I’m wondering if it could be something to do with a resampling issue like the one I mention with the R factor. If that is indeed the case, then it seems to me the same logic would apply to the per pixel values in the K factor raster, since they are also calculated as a per ha value. I could be really getting the wrong end of the stick with all of this though.
The “aligned_erodibility_[SUFFIX].tif” intermediate output is intentionally inside a directory named “/intermediate_outputs/churn_dir_not_for_humans/”, so we do not recommend attempting to interpret any files within that subfolder unless you are a machine. That being said, my understanding of SDR’s erodibility layers is that the units are in (tons * hectare * hour) / (hectare * MJ * mm), never are they given as per pixel, but someone please correct me if I’m mistaken. Further, your inclusion of MJ in the numerator leads me to wonder if you are mixing up erosivity and erodibility.
Without validating results against observed observations, we recommend interpretation of output values as relative patterns through space, and advise against putting too much confidence in the absolute numbers. Additionally, please note that the USLE equation that underlies the SDR model is known to not work very well in steep areas, perhaps such as the mountains of Lesotho. It often produces incorrectly high erosion values in such places. Please review the SDR model limitations in the User Guide to help you understand the appropriate use and interpretation of results.
I hope this makes sense. It looks like another reply is being drafted now too, which will likely be more helpful.
Thanks for your question! The calculations look correct to me. Though the pixel area changes, the resampling step should not change the erosivity values. At that step, they remain per constant unit of area (hectare). The aligned_erosivity raster units and values should be the same as the input erosivity raster.
where cell_area_ha is the pixel-to-hectare conversion factor, defined as:
cell_area_ha = cell_size**2 / 10000.0
So the units are:
LS factor * erosivity * erodibility * pixels-to-hectares conversion factor
[unitless] * [MJ * mm / (ha * hr * yr)] * [t * ha * hr / (MJ * ha * mm)] * [ha / pixel]
which cancels to t / (pixel * yr).
I would recommend spot-checking by identifying the R, K, and LS values in your intermediate outputs at a single point, then trying the math by hand. Please do let us know if this still doesn’t sound right!
You are right Jesse I had erroneously said erodibility at one point when I meant erosivity so I have edited my post accordingly.
I did look at manually calculating RKLS and it does indeed reduce the value correctly to account for the pixel size as you explained @esoth. In other words, the outputs are indeed in tons/pixel, as the user guide says. I was pretty sure this would have been integrated into the software somewhere but just got a little concerned when I looked at the aligned_erodibility.tif, which is not recommended in any case by the sound of things!
It seems like it is a case of the USLE approach over-estimating erosion in mountainous terrain. For the purposes of this relatively quick study I am doing, it seems that is just something that I must work with as an inherent limitation to the approach?
I’ll also note that the USLE equation that underlies the SDR model is known not to work very well in areas of steep slope, and this often produces incorrectly high values of erosion in these places.