Nice, this is an elegant solution. A quick Google search suggests that the most extreme slope normally found on Earth (apart from cliffs) is about 60 degrees, i.e. a cos(i) of 0.5. So any value below that is probably DEM mismatch.

I'm fine with keeping it as a config option without an explicit apply_oe.py override, since (a) the value of 0.3 seems quite conservative, and (b) anyone running apply_oe.py is already committing to a lot of our "best practice" rules of thumb.

I would like to discuss more on the cos(i) threshold .

To begin, cos(i) below 0.5 is certainly plausible and depends on sun angle. e.g., SZA=40 , SAA=180, slope=25, aspect=359 -> ~0.42.

Additionally, as a caution, this discontinuity would not allow for cos(i) to be solved for optimally from radiance. Or would just need to ensure this config is updated dynamically when solved in this mode.

Alternatively to this, I have a similar, yet slightly different idea that could be worth trying. In Richter (1998), they talk a little about this issue where for cos(i)<0.3, reflectances can "rise to a factor of 2-6". In this work, they apply a linear function for this range in reflectance space, to dampen the bright surfaces. However, my loose interpretation of this is that this could also be applied to radiance space.In that we could amplify cos(i) for these regions based on EQ-15a for just the direct radiance terms. ... In reflectance space, this is their eqn (also shown in Fig 4.):

$$ G = cos(i) / cos(B_T)$$

$$ g <= G <= 1$$

$$ \rho_g = \rho_f * G$$

In Richter 1998, they used

• $$B_T = 65 degrees$$,
• $$g$$ = 0.33

This of course is also empirical similar to yours here, @pgbrodrick but may be worth experimenting with since we have all of the available terms for this method too.

To be clear, what I am suggesting (and is NOT in Richter 1998), is to inflate cos(i) for the direct radiance terms for low incidence angles by

$$cos(i)_G = cos(i) * (1/ G) $$

My hypothesis is that this should effectively get you to a similar reflectance state, but keep the function continuous and also avoid any cutoff-artifacts that could happen from having a "threshold" based approach. Air quotes because this would also be a threshold based approach :)

@pgbrodrick Edit: contrary to everything I just said above, my statements would still be cosmetic really and just be harder to manage and turn off. ..

Your method would be much easier actually and one could solve for cosi and just set min cosi to be zero.

@brentwilder - based on your feedback and @davidraythompson's comments today, I modified the cos_i default to be 0.0, so that we're not preventing the free solve silently. Instead, I set the heuristic 0.3 threshold (subject to change) inside template_construction.

In the free-solution case, I think as long as the bound is low, the objective function should stay continuous, and we shouldn't have an issue. Given that we use a fully per-pixel cos(i), I'm not sure the above * (1/G) based manipulation is necessary, but I could be swayed! I propose this as a stopgap, that perhaps we can make better soon. Certainly closing the gap to a full cos(i) solution per-pixel is the way to go if we can get it to work universally.

Great ideas and discussion! I agree with your solution as a stopgap, @pgbrodrick. Subject to further improvements down the line. Merging for now.