Problems verifying SeaDAS implementation of QAA

[blesht]

Hi - I have converted the SeaDAS QAA code (qaa.c and get_qaa.c in l2gen msl12 9.5.0-V2019.3 (Aug 16 2019 13:00:37)) to run in R.  To be sure I got things right in the translation, I did the following experiment (based on MODIS data).  I generated a large set L2 files adding the products a_488_qaa, bbp_488_qaa, bb_488_qaa, and b_488_qaa to the output.  Then I used the Rrs_{nnn} values from that same set of files to drive my QAA.R code which outputs these same products.  In an ideal world (I thought) the values of these products from my R code would exactly (or pretty close to exactly) match those generated by SeaDAS.  Alas, that turns out not to be the case.  I've gone over my code many, many times and although there might be some bug(s) that I've missed, there are a few things about the SeaDAS output itself that I don't understand.  I'm hoping that if I describe my results someone might provide some insight into the source(s) of the discrepancies.

First:  Although the SeaDAS code get_qaa.c shows that b_488_qaa should be a linear function of bb_488_qaa, this is only true for values of bb_488_qaa < 0.107865 above which SeaDAS outputs b_488_qaa as a constant with value 5.7767.   I haven't found anything in the SeaDAS source that would account for this.  It must be there somewhere, but where (and why)?

Second:  bb_488_qaa is simply the sum of bbp_488_qaa and bbw_488.  Since bbw_488 is a constant (as noted in the reply to one of my earlier queries, https://oceancolor.gsfc.nasa.gov/forum/oceancolor/topic_show.pl?dln=28423;pid=47947) the difference (bb_488_qaa - bbp_488_qaa) should equal bbw_488.  However, it does not.  The value of bbw_488 is given (in the L2 header) as 0.001437299 but the differences range from 0.00108 to 0.00144 with the most common value (about half the observations) being 0.001435.  This may be a small point, but again, why doesn't the difference equal a single value equal to what is supposed to be the input value of bbw_488?

Finally:  When comparing my R results to the SeaDAS results, I find that they (a, bbp) agree pretty much as expected (slope~=1, intercept~=0, r^2=1) only for the cases when Rrs_667 >= 0.0015.  For cases below that Rrs_667 value the slope of my a_R vs SeaDAS a_488_qaa is 0.7948, the intercept is 0.005728 and the r^2 value is 0.957.  Beside the non-unity slope, the increased scatter seems curious insofar as both my code and SeaDAS are (presumably) using the same input Rrs spectra.   Because there is a switch in the QAA(v6) algorithm from one reference wavelength to another based on the value of Rrs_667 this was an obvious place to look for a problem in my code, but I've gone over it several time without finding an error in coding.  So, I guess I want to know if it is possible that the qaa.c source does not correspond to the qaa routine actually implemented in SeaDAS and also whether anyone else has attempted an independent verification of the SeaDAS QAA products.

Thanks!

[OB.DAAC - SeanBailey]

Barry,

1) The total scattering  product has a valid range of 0.0001 to 5.  It's stored as a scaled integer, so with the scaling applied the maximum value would be 5.7767.  My suspicion is that the pixels you're seeing pegged are probably flagged as PRODFAIL.

2) The default is to set seawater_opt=1, which enables a temperature and salinity dependent correction to bbw

3) The code listed here https://oceancolor.gsfc.nasa.gov/docs/ocssw/qaa_8c.html is the code as implemented in l2gen with the latest version of SeaDAS
The farther south of 0.0015 for Rrs_667 you get, the more influence the bbw will have on the total bb (and thus bbp), so perhaps #2 is showing it's impact here.
Try running l2gen with seawater_opt=0.

Sean

[blesht]

Ah!  Thanks (again), Sean.  I'm sure that is the explanation.  I just wish I'd noticed it before I did my own reprocessing of the entire mission (my little laptop is getting tired).  Lesson learned - read and re-read the l2gen documentation.  At least I can stop searching for a bug in my R code.

Best, Barry

[blesht]

Hi Sean - Just to add to the record here, no joy, I'm afraid.  I did a complete MODIS reprocessing with seawater_opt=0 and the results are essentially the same.  I get excellent agreement between my code and SeaDAS for the cases when Rrs_667 >= 0.0015, but when Rrs_667 < 0.0015,  my_a = 0.801*seadas_a + 0.005 with r^2=0.957 and my_bbp = 0.868*seadas_bbp - 0.00016 with r^2=0.984.   So its back to bug hunting for me.   If anything else should occur to you, please let me know, but I'm going to operate on the assumption that the problem is mine.  Thanks, Barry

[blesht]

OK - problem is solved.   It turns out that I made the incorrect assumption that the MODIS 555nm band would be the appropriate one to use for what the QAA algorithm uses as the reference wavelength for the cases when Rrs_667 < 0.0015.  In my (weak) defense, the index for the reference wavelength is referred to in the qaa.c code as "idx555" and, in the various forms of QAA documentation including published papers, 555nm is referred to as the reference wavelength for the low Rrs_667 cases.  In fact, for MODIS, SeaDAS uses the 547nm band as the lower range reference.  Had I taken a closer look at the L2 file header, I would have noticed this:      "qaa_wave =  412,  443,  488,  547,  667\n".  By using Rrs_547 rather than Rrs_555 as the reference my R computations match the QAA output exactly - whether those values are what the QAA is 'supposed' to produce is still a question, though not one I'm going to worry about now.

[OB.DAAC - SeanBailey]

Barry,
Glad you found the problem :smile:  The code does indeed refer to 555nm, so I can see that as being confusing for MODIS, which has 555nm, but it's not an "ocean" band, so the code uses 547nm.
The 555nm band is a wider band (20nm vs 10 for 547) with lower signal to noise ratio (~225 vs ~750 for the 547nm).

Sean

[blesht]

Yeah - you'd think I've been doing this long enough to have realized that in the first place.:grin:  Forging ahead now.  Thanks again (as always).   Barry