When is reflectance actually reflectance?

[djhoese]

Problem

WARNING: I'm not a scientist.

This is related to #52, specifically the standard_name used for visible channel data returned by SatPy. We return visible channel data in the "reflectance" calibration with the standard_name of toa_bidirectional_reflectance. By CF definition this term should only be used for reflectance data that includes the normalization for the solar zenith angle (SZA). This normalization (as some may know) is applied by dividing by the cos(sza). While talking with @kathys, she mentioned that she has always referred to the data without this SZA adjustment as "normalized radiances" (normalized for the theoretical maximum energy entering the atmosphere from the sun, not for the angle of the sun).

The CF definition for toa_bidirectional_reflectance:

"Bidirectional_reflectance" depends on the angles of incident and measured radiation. Reflectance is the ratio of the energy of the reflected to the incident radiation. A coordinate variable of radiation_wavelength or radiation_frequency can be used to specify the wavelength or frequency, respectively, of the radiation. "toa" means top of atmosphere. toa_bidirectional_reflectance includes a factor to account for the cosine of the solar zenith angle but does not include any integration over solid angle.

So, I'm wondering if we should be listing the pre-sunz corrected reflectances as having a different standard_name. Thoughts?

Special Notes

• VIIRS SDRs already have the sunz correction done to them so toa_bidirectional_reflectance is actually correct for that reader
• Both pre-sza and post-sza adjustment data are unitless, so this can't be used to know which is which.
• I can't find a CF standard name that seems to apply to the pre-SZA adjustment reflectance data.

[adybbroe]

@djhoese Thanks for bringing it up. First, I think the topic/name of the issue should be "When is reflectance actually a reflectance?" Right? Then, I have always thought of and when discussing with my colleagues here in Sweden understood that when we say "normalised reflectance" it is the TOA Bi-directional reflectance, thus corrected for the sun zenith. So the uncorrected we usually sloppily just call "reflectance".

Concerning finding a suitable name for the reflectance when the sun-zenith correction has been omitted, I would thus lean towards something like "unnormalised reflectance". But I would be surprised if CF has that?

[adybbroe]

In the NWCSAF/PPS software we have intermediate level-1c files storing Tbs and Reflectances in netCDF. There we do:

        description = AVHRR ch_1
        id_tag = ch_r06
        long_name = Reflectance for 0.6 micro meter channel
        scale_factor = 0.01
        standard_name = toa_bidirectional_reflectance
        sun_earth_distance_correction_applied = TRUE
        sun_earth_distance_correction_factor = 0.94129384
        sun_zenith_angle_correction_applied = TRUE
        units = percent
        valid_range = 0,20000
        wavelength_range = 0.58,0.68

So, just to be VERY clear, what we mean! :-)

[frarue]

@djhoese @adybbroe Good to thoroughly discuss such stuff (and I really like your very clear NWCSAF/PPS netcdf statements also about the sun-earth distance correction). I currently work a lot on VIS recalibration and I always refer to Nicodemus et al. p. 12: "A reflectance factor is defined as the ratio of the radiant flux actually reflected by a sample surface to that which would be reflected into the same reflected-beam geometry by an ideal (lossless) perfectly diffuse (lambertian) standard surface irradiated in exactly the same way as the sample". That said I do not understand the definition of @kathys because the "theoretical maximum energy entering the atmosphere from the sun" cannot actually be determined without considering the sun-zenith angle. Or am I wrong? Reflectance computation without SZA would be simply dividing the radiance by the solar irradiance. I cannot easily see the benefit of using this against radiance but: might be worth to consider what unit it would have: Assuming radiance comes in [Wm-2sr-1] and the solar irradiance in [Wm-2] ...it would then be [sr-1]? ...probably not a good idea to call this "unnormalized reflectance" then. My suggestion would be to either do the full reflectance computation (apart from cos(sza): don't forget π!) and call it toa_bidirectional_reflectance or to keep the radiance. But maybe I'm missing something?

[mraspaud]

Dividing the radiance by the solar irradiance has been done and used a long time in different places, so that's probably why it's still in use today. But I concur that the visible channels should be sun normalised. However, as @adybbroe noted in an offline discussion, the cos(sza) has the problem of becoming zero at the day/night terminator, so for now this is smoothed out from 88 degrees onwards in satpy iirc. And this raises then the question: is the data then truly toa_bidirectional_reflectance ? Another aspect is that we are sometimes using another normalisation accounting for a better effective solar pathlength (Article from Li and Shibata 2006). Is the resulting dataset then a toa_bidirectional_reflectance ? For imagery, we maybe don't need to be really picky on what the physical meaning of the dataset is, but for data analysis purposes it is of utmost importance. So I'm wondering now if we should rename our modifier to make things a bit clearer, and take better care of the standard name and units.

[sfinkens]

On NOAA's GOES Calibration page they call the uncorrected value "reflectance factor" or "effective albedo" with the additional explanation:

The value of 1 corresponds to the radiance of a perfectly reflecting diffuse surface illuminated at normal incidence when the sun is at its annual-average distance from the Earth.

I guess that's what @kathys wanted to express. There is no standard name for effective albedo, though.

@frarue Regarding the units: On the NOAA page the reflectance factor A is computed as A = k*R where R is the radiance in W m-2 sr-1 um-1 and k is the calibration coefficient in m2 sr um W-1, so that A is indeed unitless. I agree that this quantity is not very useful for a proper quantitative analysis. But I could imagine that for a qualitative analysis (for example by a forecaster) a value between 0 and 1 might be easier to interpret than some W m-2 sr-1 um-1 - even if the value is not perfectly accurate.

[mraspaud]

Good point about the units @sfinkens, I would even add that radiances from differents sensors/platform have often different units (some have cm instead of um from example), which doesn't make it easy to compare.

[frarue]

@sfinkens Really...they simply set SZA=0. Interesting, thank you for pointing me there! So in summary the reason why A is useful is because it normalises the radiances to the max possible for the sensor's SRF and thus get's rid of the many different units for a human observer? A sensor-normalised radiance, so to say? How about storing it as radiance with the sensor's k as attribute? That would then be CF compliant, I guess?

[sfinkens]

@frarue I'm not sure there is such a factor for each sensor. In some places I have seen lookup tables being used for calibration.

What about making both datasets available: Something like calibration='reflectance_quick' without correction and - if possible - calibration='reflectance' with solar zenith and/or pathlength correction.

[frarue]

Just one last thought that crossed my mind: The spatial pattern of the image probably is rather comparable to radiance than to reflectance (as values decrease towards pixels with higher SZA). So maybe calibration='radiance_normalised' would be more accurate than reflectance_quick?

[djhoese]

@frarue To be clear, the parts I included about my discussion with @kathys were my interpretation and how I understood them. They may not be scientifically accurate and are likely more confusing coming from me. She also likely simplified things for me to make it easier to understand for the specific discussion we were having 😉

Otherwise, here are some things that come to mind when reading this discussion:

• I forgot to link @leonhardScheck to this issue (the person who originally brought up this issue in #52)
• SatPy has, in general, been leaning towards giving the user the most flexibility when requesting the data they want to load. This may mean giving them less useful data immediately or in the default behavior, but giving them complete control over what corrections/adjustments/normalizations are applied to the data. If SatPy was only for making images this wouldn't be a big deal, but for data analysis (as @mraspaud pointed out) this is a big deal and something we would like SatPy to be a good option for.
• I think in the spirit of the above point it would be best if an SZA adjustment wasn't automatically applied to visible band data if it doesn't have to be. @mraspaud's comment about the other path-length modifier in SatPy leads me to this. If someone wants to come up with a better way of calculating this adjustment then it would be really nice if they could use SatPy to do it and make it available to everyone using SatPy.
• I don't think we can have readers only provide radiance data and not provide the "normalized radiance" (pre-sza reflectance) data since some formats have special calibration factors and/or make loading this data an in-direct process (using look-up table variables inside the file, etc).
• I think we definitely have to have a separate standard_name for the pre-SZA data. Regardless of what we think it should be, should we suggest this be added to CF? From what the rest of you have said this data isn't exactly useful for data analysis purposes?
• If we add a calibration level for this I think it should be normalized_radiance (US English, adjective first) or a similar/better name. Note that adding this may require changed the python code for a lot of readers that specifically check if ds_id.calibration == 'reflectance': # produce reflectance-like data.
• I think the conversation here and in other recent issues is making it obvious that SatPy needs a way to say "I want modifiers applied for X, Y, and Z for all products that they could apply to". This gets very complicated with compositors considering some people may not want them applied for certain variations of a compositor as these modifiers are now explicitly in the compositor's recipe. Additionally this just moves the "magic" from the readers to some where else. Even worse, the decision of what "magic" is used is in some weird place between data loading and compositor/modifier generation. This place is likely the dependency tree which is easy for me to understand (because I wrote it) but explaining it to a user and in the SatPy docs won't be so easy. It isn't technically one of the core "steps" of satpy (reading, resampling, compositing, writing).

</braindump>

[sfinkens]

A separate standard_name is probably the more pragmatic solution than adding a new calibration. Users performing a qualitative study wont't care so much. And users performing a quantitative study will probably check out the exact definition - for which the standard name is very useful.

[mraspaud]

I did some more thinking about this, reading the above comments and also #531. For me, one thing is sure is that everyone, especially scientists want to have control about what is applied to the data they work on. So instead of adding calibration levels, I would like to add modifiers. However, I understand @djhoese's concern about algorithms falling in strange locations, so to prevent that, I propose to enable the modifiers that fit inside the readers. Let me give an example with a visible channel as this is what we are talking about here: Fully calibrated and corrected channel: modifiers = (space_masked, solar_irradiance_normalized, sunz_corrected, sun_earth_distance_corrected) The get_dataset method of the file handler would get the list of desired modifiers and apply the ones it knows about, in this case space_masked and solar_irradiance_normalized. The other would be handled by "regular" modifiers as they are today.

That would require some modification in the handling of the dataset modifiers in the reader yaml files, but I don't think it'll be overly complicated.

The benefit is that almost all aspects of the dataset creation could be switched on or off at the users' request. Also, this makes the reader magic disappears :)

Another thing I'm considering is making calibration levels aliases for a given list of modifiers.

Would all of this make sense ?

[sfinkens]

Brilliant! One note: It might happen that some space pixels are masked during calibration even if the space_masked modifier is not specified. For example if a space radiance is 0 it cannot be converted to brightness temperature using the planck formula.

[mraspaud]

Yes, but that would fall into calibration-induced invalid data I suppose. It could be that some geo readers don't even need to implement space masking as it is already provided in the data file. In that case we would indicate it in the dataset description in the yaml file (like we do with sunz correction for viirs sdr data that has it already applied by default).

[sfinkens]

:+1: That sounds reasonable.

[djhoese]

Very good ideas. I'm a little worried this is going to make the dependency tree handling even more complicated and raises the level of understanding that a reader developer has to have. I could also be misunderstanding this so feel free to tell me I'm being dumb.

I'll put my questions in a numbered list so they are easier to answer:

1. What does a compositor recipe look like for an ABI natural_color with these changes? Right now it results in DatasetIDs with modifiers = ('sunz_corrected',). The dependency tree sees that it doesn't have DatasetID(..., modifiers=('sunz_corrected',)) so it looks for DatasetID(..., modifiers=tuple()). Modifiers depend on the order they are applied.
2. In the proposed system, the reader developer would have to add these modifiers to the DatasetID regardless of if they were specifically applied, right? Like does space_masked get set for VIIRS (polar-orbiter) data? Or are modifiers inconsistent between various readers even if the instruments are very similar?
3. Are readers limited to a specific set of "known" modifiers (for now the 4 you've listed)? Or can users specify any modifier?
4. What if a modifier requires another dataset to load? Or another file to be available?
5. What happens if a user does scn.load([...], modifiers=('rayleigh_corrected',))? Does that include the sunz_corrected modifier? Does it include the other modifiers that you've listed automatically? If so, how does a user turn those off?
6. Could changing the standard_name in the readers be good enough for now and in the future we implement this "default modifiers" system?
7. This sounds more and more like the "order form" concept that GeoIPS uses where a user has a file or specification of what they want when they ask for X. So if someone asks for "C01" for ABI and has provide this recipe/order form then we know that C01 should be made with the modifiers they provided there (includes enhancements, colormaps, etc).
8. Do we need to split loading in to multiple steps to make this simpler for the user to specify?
9. What if a user wants to use a specific implementation of a modifier in the reader? Does the reader's modifier system tie in to the existing modifier system? What if it needs to apply the sunz_corrected and the user wants to use the "Li and Shibata" method instead of the cos(sza) method (these are similar right?).

I could see this working better if DatasetID search didn't check for the exact ordered match of modifiers but rather doing a "contains" check for all of the requested modifiers being in the available/loadable DatasetID. The list of modifiers in the DatasetID (and therefore .attrs) would still be ordered. However, this could result in two datasets having the same set of modifiers but different final results depending on what order the modifiers were applied.

We may have to discuss this on slack to reduce the noise on this issue @mraspaud.

[mraspaud]

What I'm proposing here is a strictly ordered, complete list of modifiers, since their order matters.

1. the recipe would require DatasetID(..., modifiers('space_masked', 'solar_irradiance_normalized', 'sun_earth_distance_corrected', 'sunz_corrected')). The dependency tree will strip down the modifiers it knows about and send the rest to the reader.
2. space_masked is added to VIIRS as a convention, then others should be valid across all VIS channels/instruments
3. We limit them to known modifiers. We will go quite far with this I believe. This discussion ended up summing up what the correct definition of reflectance is, and those 4 modifiers were the ones that came up. As long as we explain clearly in the documentation that this is the definition we use, I don't see how other modifiers would be needed.
4. In this case it will be, for this dataset, an external modifier (no one in the reader)
5. scn.load([...], modifiers=('rayleigh_corrected',)) will not apply anything else than the rayleigh correction. Doesn't really makes sense, but that's what the user wants.
6. We would need to come up with a handfull of new standard names, upon which we will probably spend quite some time debating. I propose we just go ahead with the modifiers.
7. Maybe. We will have sensible defaults for when the modifiers aren't provided (None)
8. I'm not sure I understand, could you provide an example ?
9. At the moment we won't allow swapping. Sunz correction is a good example as viirs SDR data comes pre-corrected with it, but to change it for another correction, we probably have to go back to radiances and reapply calibration, etc. I think this is another issue, but worth investigation.

[mraspaud]

If we roughly agree on the concept, we could maybe open a new issue or PR to discuss the implementation details.

[gerritholl]

What is the current/previous way to get the correct reflectance where this is not built in to the reader's calibration routine? There does not appear to be a modifier sun_earth_distance_corrected at the moment, does that mean this is not currently possible? From the name effective_solar_pathlength_corrected could mean sun earth distance, although usually solar pathlength means the pathlength between TOA and ground rather than between Sun and TOA.

[pdebuyl]

Hi all,

Late to the party, but I am interested in the question. GOES and Himawari data seem already corrected. The Product User Guide for GOES states that

The Radiances product can be converted from radiances to reflectance factor or brightness temperature using information provided in the product. For the reflective bands, conversion from radiance Lν to reflectance factor ρf υ is computed as: rho f_nu = kappa L_nu where κ is the ‘kappa factor’. The kappa factor κ = ((π·d2)/E sun ) represents the incident Lambertian- equivalent radiance, where d is the instantaneous Earth-Sun distance (in Astronomical Units) ...

But reflectance should have the sunz_corrected modifier applied. Is it possible that the files are distributed already corrected (I am using EUMETCast 2km NetCDF files)? As "radiance / cos(sza)", so that the factor kappa gives the correct reflectance without having to correct for cos(sza)?

Looking at the data, it is not obvious that the data is "radiance / cos(sza)" or "radiance". In the latter case, I would expect to see a smooth radiance profile going to zero as sza increases but it does not seem to be the case.

I see that the composites are set on applying the sunz_corrected modifier for GOES in any case, so maybe the issue is settled as far as the principle of applying the modifier goes. Is it the case?