Definition and accuracy of reflectance

[LeonhardScheck]

While trying to get reflectances from MSG SEVIRI solar channel HRIT files using satpy I encountered several problems.

First of all, for calibration=='reflectance' you do not actually get reflectance according to the usual definition, but reflectance * cos(SZA) * d^-2, where SZA is the solar zenith angle and d is the sun-earth distance in astronomical units. Nevertheless the standard_name is set to "toa_bidirectional_reflectance" -- this is confusing. Moreover, reflectance * cos(SZA) * d^-2 is not really a useful quantity...

It would be great if satpy could calculate real reflectances and it would be even better if the accuracy of these reflectances was sufficiently high for applications like data assimilation.

There is some code to account for cos(SZA) in the composite module. However, this code computes SZA for a fixed time, the "start_time". Ignoring the variation of the solar zenith angle during image acquisition can cause large errors when the sun is low and causes errors of a few percent in the northern hemisphere even during the day (SEVIRI scans from south to north, so start_time is wrong by about 11 minutes in Europe). For this reason I am computing cos(SZA) like this:

mu0 = np.ma.zeros((ny,nx))
Tacq = gs.datasets[chan].info["end_time"] - gs.datasets[chan].info["start_time"]
for j in range(ny) :
    tacq = gs.datasets[chan].info["start_time"] + timedelta( seconds=(j/float(ny))*Tacq.seconds )
    mu0[j,:] = np.ma.masked_outside( cos_zen( taqu, lon2d[j,:], lat2d[j,:]),0.035, 1, copy=False)

Ignoring variations in the sun-earth distance can cause reflectance errors of up to +-3%.
To account for the variation I am now using
reflectance *= float(sunpy.sun.sunearth_distance( t=time_slot ) / astropy.units.AU)**2
Maybe sunearth_distance could be integrated in pyorbital?

Only with these two corrections I get reflectances that are compatible (error<1%) with what I got using other satellite software packages.

[adybbroe]

@LeonhardScheck Thanks for your comments and code suggestions, I think they are very relevant.
We have in general not been using satpy/mpop to get accurate reflectances in an absolute manner for quantitative applications. Normally we do correct for sun zenith to get better looking RGBs. There it is not crucial to derive the sun-zenith angle with the time valid for each individual pixel, but it would be needed for surface albedo retrievals for instance. Also, we always allow the user to do his own calibration by providing the radiances. For the same reason we have not been bothering about sun-earth distance corrections. But, I agree we could, and think we should. As, far as I remember people are not always using the same sun-earth distance correction, but as long as we thrive for the most accurate and document what is being applied it is fine I think.