Radiometric correction of rugged terrain at low sun angles

[olegzhelezn]

Hi David,

I guess it is a classic problem with most algoritms: northern slopes are poorly radiometrically corrected, whereby they either get very low BOA values, or are converted to NA. I wonder if there is any workaround for this possible in FORCE.
Left: L1, right: L2 (topo & atmo corrections).

Not using winter data is not a preferred option in this case 😉

[olegzhelezn]

Here is also, for comparison, L2 from USGS :

And NDSI calculated from L2 produced by USGS (left) and FORCE (right):

I have done some testing, and it seems like the same problem exists when I do not do any atmospheric correction in FORCE (meaning TOA). It looks equally bad:

Leads me to thinking this is some sort of rounding issue maybe?

[albaviana]

Hi Oleg, David,

We found comparable artifacts in north slopes in the Alps in winter images, running version 3.10.04:

Also, we are seeing different behavior over water bodies, with a noticeably larger amount of NoData pixels in the newer version. One difference in our setup is that we replaced the ASTER DEM with the Copernicus DEM (keeping the other parameters the same), so we are not sure whether the observed changes are version-related, DEM-related, or a combination of both.

Happy to hear what do you think! Thank you in advance!!

[davidfrantz]

Hi both,

can you please confirm whether this is indeed related to the FORCE version? I just went through the version history and the last change that could explain a different behaviour was introduced in 3,7.8 (published 05.09.2022).

In general, FORCE filters out data if the TOA reflectance is already < 0. This can happen for water and shadow targets.
The main reason for doing this, is that we usually cannot recover from this. Atmospheric correction decreases the values even more (because path reflectance etc.).

The main reason why I implemented it like this is: the topographic correction makes things worse, as it essentially is a multiplicative factor that should brighten up a dark pixel. But if our reflectance is negative, the reflectance unintuitively becomes even more negative when applying the TOC factor, hence increasing the darkening effect.

Does USGS have positive reflectance in these cases?

[olegzhelezn]

Here are example values for one of the masked by FORCE pixels:

Here are the values for a neighboring (not masked) pixel:

Last one is BOA from FORCE.
It seems from just looking at the values (tested a few pixels) like whenever one of the DN in L1 is <4, the pixel is masked out.

[davidfrantz]

I don't quite understand what which table is...

[olegzhelezn]

I can confirm different versions produce different results.
Here is 3.8.01 (left) and dev (right):

The only difference in the products is the amount of masked out pixels, otherwise everything is the same.
What is interesting is that there are also some individual dark pixels in both cases that get extremely high BOA values after correction (seen in bright blue or yellow especially visible in the left image). They are not present in TOA, so they are a result of atmosperic correction (most likely).

For comparison, the same place in USGS L2 (huge topographic shadows, but all pixels there):

[olegzhelezn]

I don't quite understand what which table is...

The names of the layers indicate what is what, first two L1 and L2 USGS, other three L1, L2 USGS and L2 FORCE.

[davidfrantz]

which image is this? And what scaling factor is USGS using for L2?

[davidfrantz]

OK got it, at least the scaling factors:

If I get this right, this also results in negative values (outside of valid range) for the USGS L2 product.

I just downloaded a Level 1 image from the Alps (November). The "usual" DNs for very dark areas are all well in the 5000s (for this image at least).

The DN to TOA conversion in Landsat images is:

ADD + MULT * DN / cos(sun-zenith)

ADD and MULT are coming from the metadata. In the November image, these are -0.1 and 2e-5, respectively, hence any DN < 5000 would become negative.

If you share the image ID, I can further investigate why there are more NAs.

[olegzhelezn]

Yes, sure, I used this one as an example: LT05_L1TP_111023_20090206_20200827_02_T1

[davidfrantz]

Got it, and at least initially, reason confirmed:

Any DN < 4 becomes negative:

# TOA reflectance of DN 4
> (-0.008007 + (2.4350E-03 * 4)) / cos((90-18.50709700)*pi/180) * 10000
[1] 54.59607
# TOA reflectance of DN 3
> (-0.008007 + (2.4350E-03 * 3)) / cos((90-18.50709700)*pi/180) * 10000
[1] -22.11566

This equation has never changed since 2016 or so. Thus, I am not sure why this should be different between versions. But I can run it through an older one tomorrow, just to confirm.

But nevertheless, like I wrote earlier, the USGS also has negative values in their product (it is just not that obvious because of their scaling). It is a winter image in the mountains, under particular poor sun illumination - plus the Landsat 5 radiometric resolution.

I can potentially allow for those negative values, but it cannot come with TOC, unfortunately.

If you would apply the TOC-factor (let's assume it would be 2) to the TOA of DN 4, you would get 109 (2 times increase, as expected).
If you would apply the same TOC-factor to the TOA of DN 3, you would get -44 (2 times decrease, surely not expected and effectively even increasing shadow effects).

I am off for today, will have a closer look tomorrow about the versions

[olegzhelezn]

Please do compare the versions, since that strangely makes a lot of difference. Thanks in advance!

Wouldn't it be possible to convert negative values to 0, for example? Or rescale them? Since there are still differences between DN 1 to 4....
Also, it seems from my sample image like the exclusion of the pixels happens always due to very low values of SWIR2. All other bands are fine, and I do not really need SWIR2 that much anyways...

[davidfrantz]

1. version tests

These are runs with 2 versions. And just to rule out some GDAL issue: once with local installation using the same (older) GDAL version for both FORCE versions, the other one using the published Docker images, where the GDAL version has changed between the FORCE versions. These runs do not use a DEM to isolate the effect of the main issue that we have negative TOA values from the beginning. These do not reproduce the issue you are describing, which makes sense since the equation and routine for converting DN-to-TOA did not change. The last column is with DEM & TOC. This time, there is a difference between the images but this is due to the CC resampling of the DEM that we introduced a while ago. It does not affect the number of nodata pixels.

This is the quantitative result of counting the number of NAs (incl. image boundary). There is no difference between versions. TOC adds a couple more NAs, but I would say in a negligible range.

run	# of NA
111023-new-docker	23783403
111023-old-docker	23783403
111023-new	23783403
111023-old	23783403
111023-old-docker-toc	23783445
111023-dev-docker-toc	23783457

Conclusion I cannot reproduce and the behaviour makes sense to me.

2. Your suggestions above

No, capping at 0 would destroy the physical coherence of the spectrum.

Rescaling doesn't solve the problem, and would introduce incorrect results when applying TOC (i.e., when TOC would be applied to rescaled reflectance with an offset). Please note that USGS is not performing a TOC at all.

3. Potential changes

The only solution that makes sense for me, would be to keep negative reflectances, but do not apply TOC if this happens. However, this would need to apply to the whole spectrum if any band would be 0 or less - otherwise the spectrum is compromised. Also note that the TOC factor is estimated in SWIR2, then propagated through the spectrum with radiative transfer.

However, I would need to re-evaluate whether keeping negative TOA values would adversely affect the internal L2 processing, or whether I would need to refactor things.

[olegzhelezn]

Hi David, this is very strange. I just rerun TOA with dev and 3.8.01 (both docker), and the results are wildly different...
The parameter file is 100% the same.

[olegzhelezn]

More specifically, in my case running 3.8.01 results in less NA than in your example above (I took a screenshot of the same area):