sertit / eoreader

Remote-sensing opensource python library reading optical and SAR sensors, loading and stacking bands, clouds, DEM and spectral indices in a sensor-agnostic way.
https://eoreader.readthedocs.io/en/latest/
Apache License 2.0
278 stars 22 forks source link

Add ST_QA band for Landsat C2L2 products #99

Closed BastienKovac closed 1 year ago

BastienKovac commented 1 year ago

When using the thermal bands from Landsat C2L2 products, the band ST_QA is a critical piece of information. This band contains the uncertainty of the Surface Temperature band (ST_B10), in K. It is used to ignore outlier values, resulting for instance from passing clouds.

There are currently, as far as I'm aware, no way to retrieve this band through EOReader so it needs to be loaded manually.

It would be interesting to be able to load this band directly. I see two ways of doing this :

  1. Add the ST_QA band to the available bands for Landsat C2L2, which is probably the most straight-forward solution.

  2. Return the ST_QA band alongside TIR_1 in the resulting XArray. The big downside of this approach is that it might break retro-compatibility, but I'd argue that the ST_QA band should pretty much always be used alongside the surface temperature.

This would also allow to pass a 'max_incertitude' kwarg to the load method, automatically replacing surface temperature pixels exceeding the given incertitude threshold by NaNs. Although this might go beyond the scope of this feature request.

BastienKovac commented 1 year ago

I just realised the TIR_1 band from EOReader does not in fact match the B10 (Surface Temperature) band, making this issue pull double duty as B10 doesn't seem to be handled either. However this band is only available for L8/9 so it would probably need a specialised subclass of LandsatProduct

remi-braun commented 1 year ago

I thanks for this issue !

I will first try to answer the second message. B10 is the band read for TIR_1, so I don't understand what makes you think the contrary 😅 Could you precise what you mean by "matching" ? Are you converting the band to BT on your side by using the Landsat CCoefficients that can be found in the metadata?

Concerning the first message: It could be a good thing to add the possibility to load masks in a not hidden way for every satellite. This would however take time. For now, you can just load this band by using the following workaround:


c2l2_prod._read_band(c2l2_prod._get_path("ST_QA"))
BastienKovac commented 1 year ago

Thanks for the quick response and the workaround, this will work wonder

About the B10 issue :

After having a closer look at EOReader's code, I got mistaken between the _map_bands_tm() and _map_bands_etm() methods, so the band is indeed correct.

The values seem slightly off to me though, which is why I initially thought of retrieving the QA band.

I first opened it and tried to convert it to Celsius by first applying B10 * 0.00341802 + 149 to go from reflectance to Kelvin, as per the product's documentation (Page 12).

If I understand correctly, the conversion between Reflectance and Kelvins is done automatically when loading the band ? However the values don't seem to match when I do the same conversion manually.

This is the result of TIR_1 - 273 (to go to Celsius) with EOReader :

eoreader

And this is the result of (B10 * 0.00341802 + 149) - 273 with Rasterio :

rasterio

Note the different values in the legends. The second plot seems more plausible to me, with it being a tile in the south of France during June 2023. It's hot down there but I sure hope we're not reaching 60°C ground temperature (yet)

remi-braun commented 1 year ago

The only process I'm doing to the BT bands are the L1 conversion (see the _to_tb function, following this tutorial)

It should be OK for L1 data, but you are right, maybe it is wrong for L2 data (to be honest, I never really use L2 data, so this may not be tested well 😓)

If I understand correctly, for L2 level, the offset and gain I should use from the MTD are the one stored in LEVEL2_SURFACE_REFLECTANCE_PARAMETERS and LEVEL2_SURFACE_TEMPERATURE_PARAMETERS, is that right ?

BastienKovac commented 1 year ago

I believe that's the case yes (I'm in no way an expert either 😅 ). The table in the documentation I linked is quite handy as a reference, here it is as an image :

image

I believe the same kind of treatment should be applied to the other bands then, but I haven't checked those. The offset and gain are also available in REFLECTANCE_MULT_BAND_{X} and REFLECTANCE_ADD_BAND_{X}

remi-braun commented 1 year ago

Yes this is what I saw. Let's correct this bug then 😓 The workaround to this is to load your band "as is" and then apply the correct coefficients:


bt = prod.load(TIR_1, to_reflectance=False) * 0.00341802 + 149
BastienKovac commented 1 year ago

Perfect ! Thanks for the reactivity

remi-braun commented 1 year ago

Could you clone the latest version and just give it a try ? I hope I haven't forgotten anything 🤞

BastienKovac commented 1 year ago

There is an issue when loading the data : The metadata can't be read.

The issue comes from the band name, which is 10 internally whereas the expected metadata name is TEMPERATURE_MULT_BANDSTB10

BastienKovac commented 1 year ago

Sorry my last comment wasn't really clear. There is a problem with the metadata not being read, which isn't caught in the try/except block because it's a TypeError, and not a ValueError :

TypeError: float() argument must be a string or a real number, not 'NoneType'

The default values are also switched, c_mul should be 0.00341802 and c_add 149.0

remi-braun commented 1 year ago

Arf sorry, I did this too fast and I haven't noticed I wasn't testing thermal bands in my CI 🙄 This should be OK now, could you please try again?

BastienKovac commented 1 year ago

Yes, all fixed. The two resulting rasters (rasterio and eoreader) are identical ! Thank you

remi-braun commented 1 year ago

Ok I'll close this issue as your demand for loading quality masks is openend in another issue: