SEVIRI day_microphysics RGB - recipe

[Miruna]

In Seviri.py, the comment on day_microphysics method says "Make a 'Day Microphysics' RGB as suggested in the MSG interpretation guide(rgbpart04.ppt)." In this guide, it looks like for Red, the data in Channel 02 (VIS 0.8) is used, but in the code of the method, I see red = np.ma.masked_where(sunzmask,self[0.8].data / costheta), 

where costheta is cos of sun zenith angle. 

Where does this costheta come from? What is the exact recipe used for Day Microphysics?

[Adam Dybbroe]

Miruna,

On 05/18/2017 03:55 PM, Miruna wrote:
> In Seviri.py, the comment on day_microphysics method says "Make a 'Day
> Microphysics' RGB as suggested in the MSG interpretation
> guide(rgbpart04.ppt)." In this guide, it looks like for Red, the data in

> Channel 02 (VIS 0.8) is used, but in the code of the method, I see /red
> = np.ma.masked_where(sunzmask,self[0.8].data / costheta)/,

> where costheta is cos of sun zenith angle.
> Where does this costheta come from? What is the exact recipe used for
> Day Microphysics?
>

The 1/cos(theta) is to convert the pseudo-reflectances in the HRIT
format to physical reflectances. The "reflectances" in the HRIT format
assumes a sun in zenith so to speak. So you need to correct for the fact
that the sun does not shine directly from above, always.

Since the derivation of 3.9 reflectance by nature gives you the true
physical reflectance (which is a physical feature of the object you are
observing) we chose to do the same for the 0.8 channel. So, all bands
have physical units.

You are right, that the interpretation guide does not mention any
sun-zenith correction, explicitly. When we implemented this RGB in mpop,
we were trying to get exact information on how it was derived in
EUMETSAT, beyond what is given in the Interpretation guide. But we could
not get a firm answer, unfortunately.

You are very welcome experimenting further, and try leave out the
sun-zenith correction. But in my opinion it gets more difficult to
interpret this RGB when you mix quantities like this. Especially at
lower sun elevations. In that case you should probably multiply the 3.9
reflectance by cos(theta).

Best regards
Adam


> --
> You received this message because you are subscribed to the Google
> Groups "pytroll" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to pytroll+u...@googlegroups.com
> <mailto:pytroll+u...@googlegroups.com>.
> For more options, visit https://groups.google.com/d/optout.


--
Adam Dybbroe,
Satellite Remote Sensing Scientist,
Numerical models and Remote Sensing,
Core Services, Swedish Meteorological and Hydrological Institute (SMHI)
www.pytroll.org
nwcsaf.smhi.se
www.smhi.se

[Miruna]

Hi Adam,

Many thanks for the quick reply and the very detailed explanation! 

Why isn't this correction also applied to the Blue channel? Is it because the 10.8 channel is infrared? I apologise if it's a silly question, my background is in computer science and I am completely new to remote sensing.

I have another question, this time on the actual computation of the solar component within the reflectance of the 3.9 channel (i.e. what we put inside the Green bean in the day_microphysics RGB). I found the computation inside the reflectance_from_tbs method in near_infrared_reflectance.py within the pyspectral package. What algorithm are you implementing?

I attach to this message some slides (also from EUM, like the rgbpart04) that I found, with the equations for calculating the solar reflectance of IR3.9. I identified in your code Equation #8 (method derive_rad39_corr), #2 and #3 (inside radiance_tb_conversion.py, the methods tb2radiance_simple and radiance2tb_simple), #6 and #7 seem to be implemented too (the denom and nomin inside reflectance_from_tbs)  but I am not sure whether these 2 exact equations are implemented...is self._solar_radiance the TOARAD from Eq 6 in the slides (The CO2-corrected, solar constant at the Top of the Atmosphere in Channel 04 (IR3.9)) or something completely different? I can't manage to find the implementation of Eq. (9) from the slides in the code. I see in the code self._solar_radiance = 1. / np.pi * self.solar_flux * mu0 so I suppose it is something different? In the code I also see no trace of the satellite zenith angle, which appears in the slides so I suppose that you are implementing another algorithm for calculating the solar component of the IR3.9 reflectance?

Thank you very much!

Best regards,

Miruna

[Adam Dybbroe]

Hi!

On 05/19/2017 05:50 PM, Miruna wrote:
> Hi Adam,
>
> Many thanks for the quick reply and the very detailed explanation!
> Why isn't this correction also applied to the Blue channel? Is it
> because the 10.8 channel is infrared? I apologise if it's a silly
> question, my background is in computer science and I am completely new
> to remote sensing.

No problem! Yes, it is because the blue band is using the 10.8
Brightness Temperature in Kelvin. And there is no reflected sunlight in
that band, it is to far from the solar spectrum.

>
> I have another question, this time on the actual computation of the
> solar component within the reflectance of the 3.9 channel (i.e. what we
> put inside the Green bean in the day_microphysics RGB). I found the
> computation inside the reflectance_from_tbs method in
> near_infrared_reflectance.py within the pyspectral package. What
> algorithm are you implementing?

Have you looked at the pyspectral documentation? There the theory behind
the derivation of the 3.9 or 3.7 (depending on the sensor) reflectance
is presented:

http://pyspectral.readthedocs.io/en/latest/37_reflectance.html

It is the same approach as the one EUMETSAT is using.

-Adam

> I attach to this message some slides (also from EUM, like the rgbpart04)
> that I found, with the equations for calculating the solar reflectance
> of IR3.9. I identified in your code Equation #8 (method
> derive_rad39_corr), #2 and #3 (inside radiance_tb_conversion.py, the
> methods tb2radiance_simple and radiance2tb_simple), #6 and #7 seem to be
> implemented too (the denom and nomin inside reflectance_from_tbs) but I
> am not sure whether these 2 exact equations are implemented...is

> *self._solar_radiance* the TOARAD from Eq 6 in the slides (The

> > an email to pytroll+u...@googlegroups.com <javascript:>
> > <mailto:pytroll+u...@googlegroups.com <javascript:>>.

> > For more options, visit https://groups.google.com/d/optout

> <https://groups.google.com/d/optout>.

>
>
> --
> Adam Dybbroe,
> Satellite Remote Sensing Scientist,
> Numerical models and Remote Sensing,
> Core Services, Swedish Meteorological and Hydrological Institute (SMHI)

> www.pytroll.org <http://www.pytroll.org>
> nwcsaf.smhi.se <http://nwcsaf.smhi.se>
> www.smhi.se <http://www.smhi.se>

[Miruna]

Hi Adam,

Many thanks for the explanation. I saw the pyspectral documentation on this topic but I think the equations are different from the ones in conversion.ppt. I am trying to implement the equations from the slides (using pytroll for loading data, reprojecting, resampling, doing band additions, generating the geotiff etc.) and compare the output with the one generated by day_microphysics but I can't manage to find in the API a way to compute the satelite zenith angle which is used in the EUM recipe for computing TOARAD (The CO2-corrected, solar constant at the Top of the Atmosphere in Channel 04  (IR3.9)). Does pyspectral provide a method for computing the satellite zenith angle?

[Adam Dybbroe]

Hi,

I am very familiar with the conversion.ppt slides. They do the same as
in pyspectral, except that we use the spectral response function and
they use the central wavelength and the to parameters A and B. Using the
relative spectral response should be more accurate/correct, and is
equally applicable to all satellites. But the difference to this more
simplified approach should really be insignificant.

Be aware that in pyspectral we do apply the co2 correction as provided
in conversion.ppt. Also, be aware that this co2 correction is very
empirical. I tried to get more information from Daniel Rosenfeld than
what is in the slides, but this is all there is it seems.

Pyorbital can give you the sun-zenith angles.

-Adam

> > www.pytroll.org <http://www.pytroll.org> <http://www.pytroll.org>
> > nwcsaf.smhi.se <http://nwcsaf.smhi.se> <http://nwcsaf.smhi.se>
> > www.smhi.se <http://www.smhi.se> <http://www.smhi.se>

> >
> > --
> > You received this message because you are subscribed to the Google
> > Groups "pytroll" group.
> > To unsubscribe from this group and stop receiving emails from it,
> send
> > an email to pytroll+u...@googlegroups.com <javascript:>
> > <mailto:pytroll+u...@googlegroups.com <javascript:>>.
> > For more options, visit https://groups.google.com/d/optout
> <https://groups.google.com/d/optout>.
>
>
> --
> Adam Dybbroe,
> Satellite Remote Sensing Scientist,
> Numerical models and Remote Sensing,
> Core Services, Swedish Meteorological and Hydrological Institute (SMHI)
> www.pytroll.org <http://www.pytroll.org>
> nwcsaf.smhi.se <http://nwcsaf.smhi.se>
> www.smhi.se <http://www.smhi.se>
>


--

[Miruna]

Dear Adam,

Thank you for the explanation. I have a very basic question regarding the format of the data initially loaded from input products. I am starting from a very basic example (according to the Pytroll quick start):

from mpop.satellites import GeostationaryFactory

from mpop.imageo.geo_image import GeoImage

import datetime

time_slot = datetime.datetime(2017,05,16,12,00)

global_data = GeostationaryFactory.create_scene("Meteosat-10", "", "seviri", time_slot)

globe=global_data.project("met09globeFull")

What physical measurement is actually inside globe[3.9].data? 

For this example, np.amax(globe[3.9].data) is  335.57027648380279

                             np.amin(globe[3.9].data is -2.9250630358043366

Seeing the code in seviri.py (with the sums and differences of channels for creating RGB composites), I thought the values inside globe['channel'].data are brightness temperatures but I think this is not it.

Best regards,

Miruna

[Adam Dybbroe]

The unit is Kelvin. The parameter is Brightness temperature.

Try:

np.histogram(globe[3.9].data, bins=100)

The mask doesn't mask out the invalid (space views) data, properly.

-Adam

On 05/29/2017 04:15 PM, Miruna wrote:
> Dear Adam,
>
> Thank you for the explanation. I have a very basic question regarding
> the format of the data initially loaded from input products. I am
> starting from a very basic example (according to the Pytroll quick start):
>

> /from mpop.satellites import GeostationaryFactory/
> /from mpop.imageo.geo_image import GeoImage/
> /import datetime/
> /time_slot = datetime.datetime(2017,05,16,12,00)/
> /global_data = GeostationaryFactory.create_scene("Meteosat-10", "",
> "seviri", time_slot)/
> /globe=global_data.project("met09globeFull")/

> Adam Dybbroe,
> Satellite Remote Sensing Scientist,
> Numerical models and Remote Sensing,
> Core Services, Swedish Meteorological and Hydrological Institute (SMHI)
> www.pytroll.org <http://www.pytroll.org>
> nwcsaf.smhi.se <http://nwcsaf.smhi.se>
> www.smhi.se <http://www.smhi.se>
>
> --

> You received this message because you are subscribed to the Google
> Groups "pytroll" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to pytroll+u...@googlegroups.com

> <mailto:pytroll+u...@googlegroups.com>.
> For more options, visit https://groups.google.com/d/optout.

[Miruna]

Dear Adam,

Thank you for the clarification.  I suppose this Brightness temperature that we load in the SatelliteInstrumentScene is the same that we would obtain if, following the conversion.ppt, I would calculate the physical radiance first (R=CAL_offset+CAL_slope*Count), then convert this to Brightness temperature (Tb=C2*vc/log(C1.....) or is it different? 

I started from these Brightness Temperature values and applied the equations from conversion.ppt for calculating the IR3.9 refl (I would like to validate the equations and to see how far off the output is from the more precise version implemented in Pytroll). 

What is the interval of values expected for the IR3.9refl? I didn't manage to find information on this. Following the equations, I obtain, for one example, values between 0.95 and 2.344 but, for the same input data, using the Pytroll day_microphysics method (and adding some prints in seviri.py to check what is obtained before inserting the data into the IR3.9refl channel), I obtain values between -5.54 and 2.68. The output image I get from the equations is very far from what it should be. Do you know if, when using the equations, I first need to apply the CO2 correction of the brightness temperature (i.e. the co2corr method in seviri.py) on the data loaded in the SatelliteInstrumentScene, prior to applying the other equations?

Many thanks for your support!

Best regards,

Miruna

[Adam Dybbroe]

Miruna,

On 05/30/2017 03:08 PM, Miruna wrote:
> Dear Adam,
>

> Thank you for the clarification. I suppose this Brightness temperature
> that we load in the SatelliteInstrumentScene is the same that we would
> obtain if, following the conversion.ppt, I would calculate the physical
> radiance first (R=CAL_offset+CAL_slope*Count), then convert this to
> Brightness temperature (Tb=C2*vc/log(C1.....) or is it different?
>

It should be the same yes. With mpop/satpy you can also choose to read
the radiances (not applying the calibration to Tbs). For instance:

glbd.load([0.8, 1.6, 3.9, 10.8, 13.2], calibrate=2)

> I started from these Brightness Temperature values and applied the
> equations from conversion.ppt for calculating the IR3.9 refl (I would
> like to validate the equations and to see how far off the output is from
> the more precise version implemented in Pytroll).
>
> What is the interval of values expected for the IR3.9refl? I didn't
> manage to find information on this. Following the equations, I obtain,
> for one example, values between 0.95 and 2.344 but, for the same input
> data, using the Pytroll day_microphysics method (and adding some prints
> in seviri.py to check what is obtained before inserting the data into
> the IR3.9refl channel), I obtain values between -5.54 and 2.68. The
> output image I get from the equations is very far from what it should
> be. Do you know if, when using the equations, I first need to apply the
> CO2 correction of the brightness temperature (i.e. the co2corr method in
> seviri.py) on the data loaded in the SatelliteInstrumentScene, prior to
> applying the other equations?
>

In theory the reflectances should be between 0 and 1 of course (or 0 and
100%). But due to 3d effects (shadows etc) they can be above 1 (100%).
Also, since the 3.9 reflectance derivation is based on assumptions that
do not hold when you have multi layer (or thin) clouds for instance you
should expect significant deviations from the theoretical case.

Thus, you can have negative reflectances. Look in pyspectral how we
filter out those unrealistic values.

I would suggest you look at histograms, rather than just min and max.

-Adam

> > Adam Dybbroe,
> > Satellite Remote Sensing Scientist,
> > Numerical models and Remote Sensing,
> > Core Services, Swedish Meteorological and Hydrological
> Institute (SMHI)
> > www.pytroll.org <http://www.pytroll.org> <http://www.pytroll.org>
> > nwcsaf.smhi.se <http://nwcsaf.smhi.se> <http://nwcsaf.smhi.se>
> > www.smhi.se <http://www.smhi.se> <http://www.smhi.se>
> >
> > --

> > You received this message because you are subscribed to the Google
> > Groups "pytroll" group.
> > To unsubscribe from this group and stop receiving emails from it,
> send
> > an email to pytroll+u...@googlegroups.com <javascript:>
> > <mailto:pytroll+u...@googlegroups.com <javascript:>>.
> > For more options, visit https://groups.google.com/d/optout

> <https://groups.google.com/d/optout>.
>
>
> --
> Adam Dybbroe,
> Satellite Remote Sensing Scientist,
> Numerical models and Remote Sensing,
> Core Services, Swedish Meteorological and Hydrological Institute (SMHI)
> www.pytroll.org <http://www.pytroll.org>

> nwcsaf.smhi.se <http://nwcsaf.smhi.se>
> www.smhi.se <http://www.smhi.se>
>


--

[miruna stoicescu]

Hi Adam,

Many thanks for all your support. I finally managed to get some nice looking output, using the equations from conversion.ppt and Pytroll for loading data and calculating sun zenith angles and satellite zenith angles.

I am using the get_observer_look to calculate the elevation and azimuth of the satellite, and I compute the satellite zenith angle using the 90-satellite elevation relation.

Regarding the arguments of this method, could you please tell me whether the altitude of the satellite should be provided in meters? 

For METEOSAT-10, I set latitude to 0, longitude to -0.4 (I found it to be 0.4 W on WMO website https://www.wmo-sat.info/oscar/satellites/view/304) and the altitude to 35786000 m according to the same website. For the lon,lat, altitude above sea level of the point on Earth, I am using the output of area.get_lonlats() and 0. Does this make sense?

Many thanks,

Miruna

You received this message because you are subscribed to a topic in the Google Groups "pytroll" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/pytroll/inmH9F4_qo0/unsubscribe.
To unsubscribe from this group and all its topics, send an email to pytroll+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Adam Dybbroe]

Miruna,

Sounds absolutely sane to me, everything!

I actually thought Met-10 was now at 0.0 degrees, but I also see the 0.4
degree, so that's probably right.

Of course one could add an elevation model and use the actual elevation
to get a more precise result instead of using just 0 as you do. But,
that seems a bit overkill for this purpose. I would have done the same
as you.

-Adam

On 06/06/2017 03:35 PM, miruna stoicescu wrote:
> Hi Adam,
>

> Many thanks for all your support. I finally managed to get some nice
> looking output, using the equations from conversion.ppt and Pytroll for
> loading data and calculating sun zenith angles and satellite zenith angles.
> I am using the get_observer_look to calculate the elevation and azimuth
> of the satellite, and I compute the satellite zenith angle using the
> 90-satellite elevation relation.
> Regarding the arguments of this method, could you please tell me whether
> the altitude of the satellite should be provided in meters?
> For METEOSAT-10, I set latitude to 0, longitude to -0.4 (I found it to
> be 0.4 W on WMO website
> https://www.wmo-sat.info/oscar/satellites/view/304) and the altitude to
> 35786000 m according to the same website. For the lon,lat, altitude
> above sea level of the point on Earth, I am using the output of
> area.get_lonlats() and 0. Does this make sense?
>
> Many thanks,
>
> Miruna
>
>
> On 31 May 2017 at 10:55, Adam Dybbroe <Adam.D...@smhi.se

> <mailto:pytroll%2Bu...@googlegroups.com>
> > <javascript:>
> > > > >
> <mailto:pytroll+u...@googlegroups.com
> <mailto:pytroll%2Bu...@googlegroups.com>

> <mailto:pytroll%2Bu...@googlegroups.com>
> <javascript:>
> > > > <mailto:pytroll+u...@googlegroups.com
> <mailto:pytroll%2Bu...@googlegroups.com> <javascript:>>.

> > Adam Dybbroe,
> > Satellite Remote Sensing Scientist,
> > Numerical models and Remote Sensing,
> > Core Services, Swedish Meteorological and Hydrological
> Institute (SMHI)
> > www.pytroll.org <http://www.pytroll.org>
> <http://www.pytroll.org> <http://www.pytroll.org>
> > nwcsaf.smhi.se <http://nwcsaf.smhi.se>
> <http://nwcsaf.smhi.se> <http://nwcsaf.smhi.se>
> > www.smhi.se <http://www.smhi.se> <http://www.smhi.se>
> <http://www.smhi.se>
> >
> > --

> > You received this message because you are subscribed to
> the Google
> > Groups "pytroll" group.
> > To unsubscribe from this group and stop receiving emails
> from it,
> send
> > an email to pytroll+u...@googlegroups.com

> <mailto:pytroll%2Bu...@googlegroups.com> <javascript:>
> > <mailto:pytroll+u...@googlegroups.com
> <mailto:pytroll%2Bu...@googlegroups.com> <javascript:>>.

> > For more options, visit
> https://groups.google.com/d/optout
> <https://groups.google.com/d/optout>
> <https://groups.google.com/d/optout

> <https://groups.google.com/d/optout>>.
>
>
> -- Adam Dybbroe,
> Satellite Remote Sensing Scientist,
> Numerical models and Remote Sensing,
> Core Services, Swedish Meteorological and Hydrological
> Institute (SMHI)
> www.pytroll.org <http://www.pytroll.org> <http://www.pytroll.org>

> nwcsaf.smhi.se <http://nwcsaf.smhi.se> <http://nwcsaf.smhi.se>
> www.smhi.se <http://www.smhi.se> <http://www.smhi.se>
>
>
>
> --

> Adam Dybbroe,
> Satellite Remote Sensing Scientist,
> Numerical models and Remote Sensing,
> Core Services, Swedish Meteorological and Hydrological Institute (SMHI)
> www.pytroll.org <http://www.pytroll.org>
> nwcsaf.smhi.se <http://nwcsaf.smhi.se>
> www.smhi.se <http://www.smhi.se>
>
> --

> You received this message because you are subscribed to a topic in
> the Google Groups "pytroll" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/pytroll/inmH9F4_qo0/unsubscribe

> <https://groups.google.com/d/topic/pytroll/inmH9F4_qo0/unsubscribe>.

> To unsubscribe from this group and all its topics, send an email to

> pytroll+u...@googlegroups.com
> <mailto:pytroll%2Bunsu...@googlegroups.com>.

> For more options, visit https://groups.google.com/d/optout
> <https://groups.google.com/d/optout>.
>
>
> --

> You received this message because you are subscribed to the Google
> Groups "pytroll" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to pytroll+u...@googlegroups.com

> <mailto:pytroll+u...@googlegroups.com>.
> For more options, visit https://groups.google.com/d/optout.

[miruna stoicescu]

Hi Adam,

Please find attached the code and an output example. You were mentioning it would be maybe useful to post the code as an example of Pytroll use. Please let me know if you see anything weird or that does not make sense in it. 

Best regards,

Miruna

 

    pytroll+unsubscribe@googlegroups.com
    <mailto:pytroll%2Bunsubscribe@googlegroups.com>.

    For more options, visit https://groups.google.com/d/optout
    <https://groups.google.com/d/optout>.


--
You received this message because you are subscribed to the Google Groups "pytroll" group.

To unsubscribe from this group and stop receiving emails from it, send an email to pytroll+unsubscribe@googlegroups.com <mailto:pytroll+unsubscribe@googlegroups.com>.

For more options, visit https://groups.google.com/d/optout.

--
Adam Dybbroe,
Satellite Remote Sensing Scientist,
Numerical models and Remote Sensing,
Core Services, Swedish Meteorological and Hydrological Institute (SMHI)
www.pytroll.org
nwcsaf.smhi.se
www.smhi.se

--
You received this message because you are subscribed to a topic in the Google Groups "pytroll" group.

To unsubscribe from this topic, visit https://groups.google.com/d/topic/pytroll/inmH9F4_qo0/unsubscribe.
To unsubscribe from this group and all its topics, send an email to pytroll+unsubscribe@googlegroups.com.