pytroll / satpy

Python package for earth-observing satellite data processing
http://satpy.readthedocs.org/en/latest/
GNU General Public License v3.0
1.07k stars 295 forks source link

When is reflectance actually reflectance? #536

Open djhoese opened 5 years ago

djhoese commented 5 years ago

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

adybbroe commented 5 years ago

@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 commented 5 years ago

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 commented 5 years ago

@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 commented 5 years ago

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 commented 5 years ago

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 commented 5 years ago

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 commented 5 years ago

@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 commented 5 years ago

@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 commented 5 years ago

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 commented 5 years ago

@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:

</braindump>

sfinkens commented 5 years ago

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 commented 5 years ago

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 commented 5 years ago

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 commented 5 years ago

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 commented 5 years ago

:+1: That sounds reasonable.

djhoese commented 5 years ago

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 commented 4 years ago

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 commented 4 years ago

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

gerritholl commented 4 years ago

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 commented 2 years ago

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?

gerritholl commented 10 months ago

Possibly related reading:

Schaepman-Strub, G., M. E. Schaepman, T. H. Painter, S. Dangel, and J. V. Martonchik. 2006. "Reflectance Quantities in Optical Remote Sensing-Definitions and Case Studies.", Remote Sensing of Environment 103 (1): 27-42. doi:10.1016/j.rse.2006.03.002, https://www.sciencedirect.com/science/article/pii/S0034425706001167

pdebuyl commented 4 weeks ago

Since my last comment, we (at RMIB) added the division by cos(SZA) to the data in our internal codebase. Empirically it seems necessary :-) It is also used in satpy: The natural color composite for ABI includes the sunz_corrected modifier for the three bands (C05, C03, C02).

The PUG defines the reflectance "factor rho f_nu = kappa_0 L_nu" where "rho f_nu" is a single quantity (not a product). L_nu is the radiance, kappa_0 = pi d^2 / E_sun. d^2 accounts for the distance to the sun, E_sun (band-dependent) for the incoming energy.

I believe that "rho f_nu" is "simply" "the reflectance that you measure from the p.o.v. of the satellite", so needing the correction.

I only found one explicit reference to this in in this document: https://www.star.nesdis.noaa.gov/goesr/docs/ATBD/Imagery.pdf

There, Eq. 3-2 confirms that rho_nu = rho f_nu / cos(SZA)

The only really annoying thing is that rho f_nu is not formally defined in the GOES PUG.

In any case, regarding satpy: ABI band data is the "reflectance factor". The calibration done in the reader reproduces the kappa_0 factor (actually, it could directly read kappa_0 in the file but recompute pi * d^2 / E_sun) and requires the cos(sza) correction for being a "toa_bidirectional_reflectance".

gerritholl commented 4 weeks ago

Just a reminder that satpy has two classes that apply a solar zenith correction:

As Schaepman-Strub et al. (2006) point out: toa_bidirectional_reflectance is not a measurable quantity, but can only be estimated using assumptions. The measured quantity is Hemispherical–conical reflectance. One of the consequences of measuring one but assuming the other is that the value can be larger than one, which is not physically possible for a real reflectance.

mraspaud commented 4 weeks ago

I think we reached a consensus yesterday at the pytroll community meeting that without the cos(sun_zen) normalisation, the quantity is not a toa_bidirectional_reflectance.

  1. So we need do find a name for that quantity, and possibly propose it as a standard name in the cf conventions. "Normalised radiance" was proposed, or maybe something along the line of "TOA radiance normalised with solar flux/irradiance" ?
  2. Assuming we can all agree that toa_bidirectional_reflectance is a good enough name for the cos(sun_zen) normalised quantity, we should make sure that requesting reflectance in satpy does in fact apply the normalisation, and make sure the other quantity is also available as another "calibration level". I know that it's probably not the exact right name, and that this quantity is probably not very physical, but it's what we can reasonably work with in the context of satpy.
gerritholl commented 4 weeks ago

I agree that we reached this consensus.

Point 1, naming is hard. How about zenith_equivalent_reflectance? What we are reporting is the reflectance as if the Sun were at the zenith (I don't remember if the actual solar irradiance is assumed or measured).

On point 2, for pragmatic reasons, I agree that we should not try to be holier than the pope. We should call it toa_bidirectional_reflectance.

Thirdly, we should document this somewhere, with a note on how we define toa_bidirectional_reflectance (here we can explain why it can be larger than 1, because it's not strictly the reflectance).

pdebuyl commented 4 weeks ago

Regarding 2, I support the idea of keeping toa_bidirectional_reflectance. In practice, any data file containing that quantity should have a reference explaining how it was obtained.

simonrp84 commented 3 weeks ago

We need to be quite careful on this, as there will be users (especially science users doing retrievals) who expect reflectances not to have the cos(sza) scaling. If this discussion results in a wholesale change to what the reflectance calibration term means in satpy then this will have to be communicated very clearly!

Regarding the naming, toa_bidirectional_reflectance seems fine for the normalised term. It's not strictly correct but it's what everyone uses. For the current "reflectance" calibration in satpy it could be renamed to toa_reflectance_factor, which I think is what GOES calls it.

mraspaud commented 3 weeks ago

I agree with the communication part, my idea at the moment is to slowly deprecate the reflectance calibration when we update satpy for this and introduce a clearer toa_bidirectional_reflectance calibration instead, along with eg toa_reflectance_factor for the non-normalized calibration.

pdebuyl commented 3 weeks ago

I like "reflectance factor" as it it what GOES uses. But is it a canonical (or consensual / obvious) way to define it?

gerritholl commented 3 weeks ago

@simonrp84 What users expect that reflectance does not have the cos(sza) scaling but toa_bidirectional_reflectance does? Reflectance is ambiguous in any case, but without dividing by cos(sza), the term reflectance is unambigiuously incorrect, so science users should not expect this to be absent for anything called reflectance. If users don't want radiances to be divided by cos(sza) they should load measurements calibrated as radiance, not reflectance.

simonrp84 commented 3 weeks ago

What users expect that reflectance does not have the cos(sza) scaling

Users who are already using satpy ;-)

Otherwise, I in the past have used satpy's reflectance (without the cos(sza) ) for BRDF modelling, simply because it's a handy way to account for solar irradiance without having to do it myself but I don't want to apply the cos(sza) as that's part of the BRDF modelling procedure anyway.

The same within some of the cloud retrieval techniques I've seen - they use satpy to load "reflectance" but correct for solar zenith themselves at another stage of their algorithm. If they had that term already applied in the data loaded via satpy, they'd have to "unapply" it.

Overall, I don't think it's up to us as developers to impose these things on users simply because it's not the 100% scientifically correct definition. We can't forsee all use cases and we don't know what shortcuts users have taken for simplicity, speed of processing, etc.

gerritholl commented 3 weeks ago

Based on Schaepman-Strub et al., I do not believe that "reflectance factor" is correct for what we are reporting:

The reflectance factor is the ratio of the radiant flux reflected by a sample surface to the radiant flux reflected into the identical beam geometry by an ideal (lossless) and diffuse (Lambertian) standard surface, irradiated under the same conditions as the sample surface.

Simon wrote:

Overall, I don't think it's up to us as developers to impose these things on users simply because it's not the 100% scientifically correct definition. We can't forsee all use cases and we don't know what shortcuts users have taken for simplicity, speed of processing, etc.

Certainly, we cannot foresee all use cases. Fortunately, Satpy already allows for all use cases by offering users (typically) a choice of three quantities when reading radiances in the solar spectrum (counts, radiance, reflectance).

But if we return a quantity using a standard term, it should ideally match the established definition for that standard term. If we return a quantity that is not described by any standard term, we should pick a term that is not yet taken and possibly propose it as a new standard term. If there is a widely used but incorrect term in the community, we can debate if we prefer to be correct or follow the crowd. AFAIK what Satpy currently returns is neither widely used, nor is there a standard term for it.

gerritholl commented 3 weeks ago

As I see it, options for the way forward of handling the default return quantity:

  1. Status quo. Continue returning a quantity called reflectance that includes the solar flux but no SZA correction.
  2. Change what we return. This would strongly break backwards compatibility.
  3. Do not change what we return, but change what we call it. This would mildly break backwards compatibility.

If we choose option (0) or (2), we can add a new calibration that returns a reflectance including sunz_correction applied.

If we choose option (1), we should add the status quo under a new name.

No matter what we call it, the status quo should remain an option.

I would personally prefer option (2). I think option (1) is too disruptive, and option (0) will continue to lead to confusion.

ameraner commented 3 weeks ago

I would also be in favour of 2, and the name change will likely have relatively minor side effects. My biggest motivation to avoid 1, apart from backwards compatibility, is that the computation of the cos(sza) correction requires an extra step with a per-pixel multiplication factor, which increases the computational effort quite a bit. I also support adding the corrected quantity as a new calibration level, since this is quite straightfoward/user friendly to request in a load call, while using the modifier is less so.

djhoese commented 3 weeks ago

I'm not sure if Satpy will support this out of the box, but I propose the below which is a detailed version of Gerrit's option 2 above:

  1. Add a new calibration level to the default set of calibrations that properly describes what we are calling "reflectance" currently. So this would be after the radiance -> reflectance calibration/normalization, but before any / cos(SZA). Examples would be normalized_radiance or reflectance_factor or something as long as it is accurate and decently well used/known.
  2. Slowly move reader by reader to use this new calibration as well as use a different CF standard_name.
  3. Update every reader to map a request for calibration="reflectance" to calibration="<new name>" and emit a warning.
  4. Update the (currently) 2 reflectance normaliation modifiers (sunz corrected and the pathlength one) to update calibration in the return DataArray attrs to be reflectance and standard_name to be toa_bidirectional_reflectance. The _satpy_id may also need to be updated, but I think this happens automatically with how modifiers work already.
  5. Release the above before Satpy 1.0 and then in Satpy 1.0 remove the deprecation and DataQuery reflectance -> normalized_radiance magic (described in 3 above).
  6. Save any Scene.load([...], calibration="reflectance") automatically including / cos(SZA) for Satpy 2.0 or a later 1.x release.
simonrp84 commented 3 weeks ago

FWIW, I asked some colleagues who know about this type of thing and they suggested that toa_bidirectional_reflectance and toa_bidirectional_reflectance_factor are suitable names for the SZA corrected and non-corrected data respectively.

gerritholl commented 2 weeks ago
  1. Update every reader to map a request for calibration="reflectance" to calibration="<new name>" and emit a warning.

This warning should also be emitted if no calibration is requested at all, at least until point 6 is implemented. Most users will directly load .load(["vis_06"]) or other solar channels, and might believe this is a reflectance, even if it's called "reflectance_factor" or "toa_bidirectional_reflectance_factor".

FWIW, I asked some colleagues who know about this type of thing and they suggested that toa_bidirectional_reflectance and toa_bidirectional_reflectance_factor are suitable names for the SZA corrected and non-corrected data respectively.

I am not convinced. The term "bidirectional reflectance factor" is widely used, and estimating this would include SZA correction. I think the term "factor" is too occupied to use for the term we are using.