Bayesian atmospheric correction over land: Sentinel-2/MSI and Landsat 8/OLI

Yin, Feng; Lewis, Philip E.; Gómez-Dans, Jose L.

doi:https://doi.org/10.5194/gmd-15-7933-2022

Articles | Volume 15, issue 21

https://doi.org/10.5194/gmd-15-7933-2022

© Author(s) 2022. This work is distributed under
the Creative Commons Attribution 4.0 License.

https://doi.org/10.5194/gmd-15-7933-2022

© Author(s) 2022. This work is distributed under
the Creative Commons Attribution 4.0 License.

Articles | Volume 15, issue 21

Development and technical paper

|

07 Nov 2022

Development and technical paper |

| 07 Nov 2022

Bayesian atmospheric correction over land: Sentinel-2/MSI and Landsat 8/OLI

Feng Yin, Philip E. Lewis, and Jose L. Gómez-Dans

Download

Final revised paper (published on 07 Nov 2022)
Preprint (discussion started on 02 May 2022)

Interactive discussion

Status: closed

RC1:
'Comment on egusphere-2022-170', Anonymous Referee #1, 09 Jun 2022

General comments:

Atmospheric correction is an essential part of satellite remote sensing of land surface. Yin et al. describe and evaluate the Sensor Invariant Atmospheric Correction (SIAC) algorithm for atmospheric correction. The existing atmospheric correction methods can be improved and therefore further development of the algorithms is welcome.

SIAC is capable of carrying out atmospheric correction both for Sentinel 2 and Landsat 8 satellite data. Furthermore, SIAC is a Bayesian (statistical) algorithm so it can take advantage of prior information and it produces uncertainty estimates for the surface reflectances produced. The algorithm is tested and validated with ground-based AERONET and RadCalNet data. There are no significant steps taken in new method development in this work, but SIAC combines well existing methods. The results shown in the manuscript show that SIAC is capable of carrying out good quality atmospheric correction.

The selections made in the development of SIAC are mostly well justified and based on previously published literature. It is very good that the authors have shared the codes for others to be used.

My main criticism is in the presentation. The quality of presentation in the manuscript varies. Time to time the text is well written and smooth but quite often the text is difficult to read. The manuscript is quite long and structured so that first SIAC and results are explained in general terms followed by the Discussion and Conclusions, and then all the technical details are mostly included in the Appendices. As this is a method development manuscript, I find this a bit difficult for the reader as it is needed to browse back and forth while reading. Furthermore, the manuscript heavily relies on use of acronyms and symbols. It is good that the main symbols are explained in a table but this also makes the manuscript very slow to read, especially for a person who is not that familiar with the field. I would strongly recommend the authors to think the use of acronyms (even single letter acronyms) and symbols, and possibly shorten and re-organise the manuscript for improved readability. Also, there were some typos in the text so proofreading is recommended. Most of the font sizes in figures are too small and very difficult to read.

Specific comments:

l.15 Abstract mentions efficient emulators. However, in the text this was a bit unclear how these were used. Could it be clarified?

l.53 AERONET is based on remote sensing, not in-situ measurements.

l.83 ".., we calculate:". "Calculate atmospheric parameter estimation"?

l.84 The list is difficult to follow.

l.129 In the abstract this was probably mentioned as "statistical emulation", neural networks are not really statistical emulation.

l.133 "This may cause errors". Can you estimate the significance of these errors

l.169 "Matrix D" was this defined?

l.183 "We assume that mean atmospheric parameters...are constant..." Can you estimate the significance of this assumption?

l.221 Can you give an example of "other artifacts"

l.223 Why the tolerance of 10% was selected? How this filtering affects the evaluation of the results?

l.300 "corrected to R", what is R?

l.333 Why different TCWV gamma values are used for S2 and L8?

p.23 Figure 13. What are the colors of the bars?

l.374 Is IQR defined?

l.375 Validation of surface reflectance uncertainty was a bit unclear. Should be clarified.

l.390 In many sentences (especially on page 25) you use "that" & "this" and it is unclear to which word these are referring to

p.26 Figure 15. What are the colors?

Citation: https://doi.org/10.5194/egusphere-2022-170-RC1
- AC2: 'Reply on RC1', Feng Yin, 05 Jul 2022
  
  We thank the reviewer for the constructive and insightfull commments. We have added our detailed reponse in the supplement.
  
  Citation: https://doi.org/10.5194/egusphere-2022-170-AC2
CEC1:
'Comment on egusphere-2022-170', Juan Antonio Añel, 15 Jun 2022

Dear authors,
After checking your manuscript, it has come to our attention that it does not comply with our Code and Data Policy.

https://www.geoscientific-model-development.net/policies/code_and_data_policy.html

You have archived your code on GitHub. However, GitHub is not a suitable repository. GitHub itself instructs authors to use other alternatives for long-term archival and publishing, such as Zenodo. Therefore, please, publish your code in one of the appropriate repositories, and include the relevant primary input/output data. In this way, you must include in a potential reviewed version of your manuscript the modified 'Code and Data Availability' section, the DOI of the code (and another DOI for the dataset if necessary).
Please, reply as soon as possible to this comment with the new links and DOIs, and be aware that failing to comply with this request will result in the rejection of your manuscript.
Juan A. Añel

Geosci. Model Dev. Exec. Editor

Citation: https://doi.org/10.5194/egusphere-2022-170-CEC1
- AC1: 'Reply on CEC1', Feng Yin, 16 Jun 2022
  
  Dear Juan,
  Thanks for pointing out the issue. I have uploaded the code into Zenodo along with the validaiton dataset used in this manuscript. The DOI to the code: https://doi.org/10.5281/zenodo.6651964 and the validation datasets: https://doi.org/10.5281/zenodo.6652892.
  Best regards,
  Feng.
  
  Citation: https://doi.org/10.5194/egusphere-2022-170-AC1
RC2:
'Comment on egusphere-2022-170', Hankui Zhang, 24 Jun 2022
I reviewed the manuscript years ago submitted to another journal + I did look forwarding to the revision which I did not see. Now I saw it in a different journal with significantly improvement in the experimental design and presentation style.

I recommended for publication as the study

presented a novel way for medium resolution (Landsat-8 and Sentinel-2) satellite data atmospheric correction

showed the presented method significantly improves the surface reflectance and aerosol property retrieval accuracy

systematically validated the presented method using Landsat and Sentinel-2 data over globally distributed AERONET

Impressive work and I have to say this is a hard paper to review as there are many techniques used in the study. I appreciated the authors make the codes public available. I like many ideas in the paper but in particular the spatial smoothness prior and the validation of the uncertainty (many geostationary satellite Bayesian AOD estimation output uncertainty without validation). I really hope the authors later (not in this manuscript) can experimentally show that the proposed method is better than LaSRC and Sen2Cor – by running LaSRC and Sen2Cor software on the datasets used in this study (this maybe boring to the authors but I believe it is beneficial to the community and to the authors if the algorithm more convinced). Or simply participate ACIX and hope to see the algorithm will be an unambiguous winner in ACIX.

I have a few comments before the publication of the manuscript

Major issues

1, The readers may wonder what is the computational cost of the proposed method and the iteration nature of the gradient descent (L-BFGS-B algorithm) may boost the computation time rapidly. Add a paragraph for discussion.

2, The authors omitted the deep blue bands for both Landsat-8 and Sentinel-2 and the band is most sensitive to aerosol among all the bands, any particular reason for doing so. MODIS does not have this band? use MODIS blue to replace not work ?

3, Spatial smoothness is interesting but unclear mostly related to how D matrix is derived. Is the spatial difference applied to the intermediate AOD and WV results derived for each iteration of L-BFGS-B algorithm. If so, Cxb is changing with iterations? The spatial difference is calculated over neighbor 8 pixels or over an entire 40 km window? How the row and column direction are averaged? Simply give an example of deriving D would solve my concern. After the line of 775, “over the whole of Gc” or of sub-40 km? Looks like the sub-40 km is used as a block to solve the equation which means the solution of each sub-40 km is independent to the neighbor sub-40 km block (but I am not sure).

The study need a consistent check of the writing style. For example, lines 102-103.

Minor issues

Discuss the possible shift from MODIS to VIIRS considering the dying of MODIS.

Line 105, how to treat b09. Is MODIS band 2 reflectance (after D2 transformation) used as a prior for b09 reflectance?

Line 675, there is another MCD43 PSF estimation paper worth citing (Che et al. 2021).

Che, X., Zhang, H. K., & Liu, J. (2021). Making Landsat 5, 7 and 8 reflectance consistent using MODIS nadir-BRDF adjusted reflectance as reference. Remote Sensing of Environment, 262, 112517.

Line 130, not quite familiar with the L-BFGS-B algorithm. But does the L-BFGS-B algorithm have its own way to calculate the gradient or the authors really used the back propagation of ANN to calculate ANN (ANN uses an automatic differentiation method) to feed into L-BFGS-B algorithm. Or back prorogation gradient is only used in B4 for RBF uncertainty but not used in L-BFGS-B?

Line 140, the LaSRC cannot be considered as DDV method anymore as each single pixel on the global have a visible/SWIR ratio (not just DDV pixels) see Vermote et al. 2016.

Line 155, the framework can tolerate incomplete coverage of Y^, how? If incomplete Y^, then the incomplete inversion of aod and wv will be spatially interpolated? Or get value using spatial smooth prior.

Line 200, “a estimates” – grammar

Line 260, how the σxaero is set

Line 310, which TOA band, the blue band should definitely greater than 10% mostly

Line 330 we are already in Section 4.3

Please label and re-arrange figures 10-12 correctly. Some are clearly surface but labelled as TOA. Check the labels of the plots. Rearrange them so that the TOA spectra and TOA plot are at the same row.

Line 470, yes, in many Landsat calibration literature, the calibration accuracy is 3%. Should this linked to appendix H?

Line 580. There is no Schaff 2021 in reference list
Citation: https://doi.org/10.5194/egusphere-2022-170-RC2
- AC3: 'Reply on RC2', Feng Yin, 05 Jul 2022
  
  We thank the reviewer for the constructive and insightfull commments. We have added our detailed reponse in the supplement.
  
  Citation: https://doi.org/10.5194/egusphere-2022-170-AC3

Peer review completion

AR: Author's response | RR: Referee report | ED: Editor decision | EF: Editorial file upload

AR by Feng Yin on behalf of the Authors (05 Jul 2022) Author's response

EF by Polina Shvedko (06 Jul 2022) Manuscript Author's tracked changes

ED: Referee Nomination & Report Request started (07 Jul 2022) by Le Yu

RR by Hankui Zhang (17 Jul 2022)

Suggestions for revision or reasons for rejection

Thanks for the authors’ response and addressing for most of my concerns. I still could not follow how the spatial smoothness prior works.

I am curious what is the spatial support in solving the Bayesian Equation (1). If there is no spatial smoothness, the Bayesian Equation can be solved at each Landsat/Setinel-2 pixel as neighbor Landsat/Setinel-2 pixels are independent. However, with spatial smoothness the Bayesian equation cannot be solved for each single pixel anymore, and a cluster of pixels will be put tougher so that D can be applied (I assume D matrix is applied on the water vapor and aerosol images). It looks like the authors choose to solve all Landsat/S2 pixels within a MODIS pixel together, i.e., D is applied all L8/S2 pixels within a MODIS pixel. In another word, the solution induced by the spatial smoothness is only applied within a MODIS pixel. The L8/S2 pixels coming from neighbor MODIS pixels are not constrained by the spatial smoothness even the two L8/S2 may be close to each other at the MODIS pixel boundaries (of course the neighbor MODIS pixels may very likely use the same prior as the prior is defined at 40 km resolution).

The above is only my guess (pity that I can only guess), but please clarify further to make it clearer and maybe pay attention to a few sentences that may cause misleading:

Line 80: “with an associated sub-40 km spatial correlation constraint” – But if my above statement is correct the spatial correlation is only on 500 m.

Table 2, Should X and Xb defined for L8/S2 resolution rather than MODIS resolution.

After Equation B3, “where n is the number of rows (columns) in the MODIS spatial grid Gc as the term is applied in the row (column) direction”. This is where my interpretation coming from since n is the number of L8/S2 pixels in a MODIS grid Gc. It implies the D with n*n size is applied to all L8/S2 pixels in a MODIS grid.

Line 855, “expressed in 6 and controlled by allowing for gradual changes in atmospheric parameters in X over the whole of Gc”. This is also where my interpretation coming from since the parameters X over Gc has gradual changes and consequently X should be defined at L8/S2 resolution (Gm).

If I am wrong, please find a place to clarify, better in 2.2 and 2.4 both.

Hide

RR by Anonymous Referee #3 (19 Aug 2022)

Suggestions for revision or reasons for rejection

The manuscript of Feng Yin proposes a method for atmospheric correction for radiometers with a typical spatial resolution of few tenths of meters relying on prior information determined at much lower spatial resolution (few hundredths of meters). The proposed method relies on a Bayesian approach. Though the title of the manuscript “BAYESIAN ATMOSPHERIC CORRECTION OVER LAND: SENTINEL-2/MSI AND LANDSAT 8/OLI” is very appealing, this research work misses delivering the expected Analysis Ready Data outcome, essentially due to several unjustified assumptions and approximations. As stated in the abstract, mitigating the impact of atmospheric effects on optical remote sensing data is critical for monitoring intrinsic land processes. In the context of delivering Analysis Ready Data, intrinsic land processes mean that no additional variables are needed to analyse the delivered surface reflectance. Despite the apparent complexity of the proposed SIAC Bayesian approach, the delivered geophysical variable, namely the bottom of atmosphere Bidirectional Reflectance Factor (factor and not function as stated in line 29) is erroneously presented as Analysis Ready Data. It still depends on the actual state of the atmosphere!

In line 60, it is written “We wish to estimate the probability distribution function (PDF) of BOA spectral BRF, R with illumination and viewing vectors Ωs , Ωv respectively”. The Bidirectional Reflectance Factor is indeed an intrinsic property of the surface that depends on Ωs , Ωv. The BRF should be estimated at the surface to achieve this objective, removing all possible atmospheric effects as written in line 94. In line 120 Equation (3), it is however assumed that the surface reflectance rb represents a “Lambertian reflectance” which is equivalent to assuming that the sky radiation is perfectly diffuse which is true only when the sky is overcast. rb is however estimated with Equation (4) which provides the BRF for the Ωs , Ωv geometry. It is therefore not a Lambertian Equivalent Reflectance as assumed in Equation 3. The relationship between rb expressed in Equations (3) and (4) is thus inconsistent. Assuming that rb in Equation (3) can be estimated with Equation (4) is strictly correct for the direct contribution only. It is inconsistent with the definition of the terms representing the total (direct and diffuse) downwelling and upwelling atmospheric transmittance. Consequently, the retrieved surface reflectance still depends on the actual atmospheric state (diffuse contribution). Hence, as long as this issue is not properly solved, the proposed approach will not work for the delivery of Analysis Ready Data. It is a major flaw in the proposed approach that should be corrected, and they are elegant ways to do so! Confusing symbols have been used to hide this issue:
- In lines 19 and 20, R defines a correlation coefficient.
- In Table 1, the symbol r is not defined.
- In Table 2, R is defined as R~N(r,Cr) as the a posteriori PDF of the BRF but the symbols r and Cr are not defined
- In Table 2, Rb is defined as Rb~(rb,Cb) as the a priori PDF of the BRF but the symbols rb and Cb are not defined
- In line 200, r is defined as the mean surface reflectance.

Sections 3 and 5 have not been reviewed as the proposed approach is erroneous. Unfortunately, I would not recommend the publication of this manuscript as long as it relies on erroneous assumptions.

Hide

ED: Reconsider after major revisions (19 Aug 2022) by Le Yu

AR by Feng Yin on behalf of the Authors (14 Sep 2022) Author's response Author's tracked changes Manuscript

ED: Referee Nomination & Report Request started (15 Sep 2022) by Le Yu

RR by Anonymous Referee #2 (03 Oct 2022)

ED: Publish subject to technical corrections (03 Oct 2022) by Le Yu

AR by Feng Yin on behalf of the Authors (10 Oct 2022) Author's response Manuscript

Short summary

The proposed SIAC atmospheric correction method provides consistent surface reflectance estimations from medium spatial-resolution satellites (Sentinel 2 and Landsat 8) with per-pixel uncertainty information. The outputs from SIAC have been validated against a wide range of ground measurements, and it shows that SIAC can provide accurate estimations of both surface reflectance and atmospheric parameters, with meaningful uncertainty information.