Source code for satpy.modifiers.atmosphere

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright (c) 2020 Satpy developers
#
# This file is part of satpy.
#
# satpy is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version.
#
# satpy is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
# A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along with
# satpy.  If not, see <http://www.gnu.org/licenses/>.
"""Modifiers related to atmospheric corrections or adjustments."""

import logging

import dask.array as da
import numpy as np
import xarray as xr

from satpy.modifiers import ModifierBase
from satpy.modifiers._crefl import ReflectanceCorrector  # noqa
from satpy.modifiers.angles import compute_relative_azimuth, get_angles, get_satellite_zenith_angle

logger = logging.getLogger(__name__)


[docs] class PSPRayleighReflectance(ModifierBase): """Pyspectral-based rayleigh corrector for visible channels. It is possible to use ``reduce_lim_low``, ``reduce_lim_high`` and ``reduce_strength`` together to reduce rayleigh correction at high solar zenith angle and make the image transition from rayleigh-corrected to partially/none rayleigh-corrected at day/night edge, therefore producing a more natural look, which could be especially helpful for geostationary satellites. This reduction starts at solar zenith angle of ``reduce_lim_low``, and ends in ``reduce_lim_high``. It's linearly scaled between these two angles. The ``reduce_strength`` controls the amount of the reduction. When the solar zenith angle reaches ``reduce_lim_high``, the rayleigh correction will remain ``(1 - reduce_strength)`` of its initial reduce_strength at ``reduce_lim_high``. To use this function in a YAML configuration file: .. code-block:: yaml rayleigh_corrected_reduced: modifier: !!python/name:satpy.modifiers.PSPRayleighReflectance atmosphere: us-standard aerosol_type: rayleigh_only reduce_lim_low: 70 reduce_lim_high: 95 reduce_strength: 0.6 prerequisites: - name: B03 modifiers: [sunz_corrected] optional_prerequisites: - satellite_azimuth_angle - satellite_zenith_angle - solar_azimuth_angle - solar_zenith_angle In the case above, rayleigh correction is reduced gradually starting at solar zenith angle 70°. When reaching 95°, the correction will only remain 40% its initial strength at 95°. """ def __call__(self, projectables, optional_datasets=None, **info): """Get the corrected reflectance when removing Rayleigh scattering. Uses pyspectral. """ from pyspectral.rayleigh import Rayleigh projectables = projectables + (optional_datasets or []) if len(projectables) != 6: vis, red = self.match_data_arrays(projectables) sata, satz, suna, sunz = get_angles(vis) else: vis, red, sata, satz, suna, sunz = self.match_data_arrays(projectables) # First make sure the two azimuth angles are in the range 0-360: sata = sata % 360. suna = suna % 360. # get the dask array underneath sata = sata.data satz = satz.data suna = suna.data sunz = sunz.data ssadiff = compute_relative_azimuth(sata, suna) del sata, suna atmosphere = self.attrs.get("atmosphere", "us-standard") aerosol_type = self.attrs.get("aerosol_type", "marine_clean_aerosol") reduce_lim_low = abs(self.attrs.get("reduce_lim_low", 70)) reduce_lim_high = abs(self.attrs.get("reduce_lim_high", 105)) reduce_strength = np.clip(self.attrs.get("reduce_strength", 0), 0, 1) logger.info("Removing Rayleigh scattering with atmosphere '%s' and " "aerosol type '%s' for '%s'", atmosphere, aerosol_type, vis.attrs["name"]) corrector = Rayleigh(vis.attrs["platform_name"], vis.attrs["sensor"], atmosphere=atmosphere, aerosol_type=aerosol_type) try: refl_cor_band = corrector.get_reflectance(sunz, satz, ssadiff, vis.attrs["name"], red.data) except (KeyError, IOError): logger.warning("Could not get the reflectance correction using band name: %s", vis.attrs["name"]) logger.warning("Will try use the wavelength, however, this may be ambiguous!") refl_cor_band = corrector.get_reflectance(sunz, satz, ssadiff, vis.attrs["wavelength"][1], red.data) if reduce_strength > 0: if reduce_lim_low > reduce_lim_high: reduce_lim_low = reduce_lim_high refl_cor_band = corrector.reduce_rayleigh_highzenith(sunz, refl_cor_band, reduce_lim_low, reduce_lim_high, reduce_strength) proj = vis - refl_cor_band proj.attrs = vis.attrs self.apply_modifier_info(vis, proj) return proj
[docs] def _call_mapped_correction(satz, band_data, corrector, band_name): # need to convert to masked array orig_dtype = band_data.dtype band_data = np.ma.masked_where(np.isnan(band_data), band_data) res = corrector.get_correction(satz, band_name, band_data) return res.filled(np.nan).astype(orig_dtype, copy=False)
[docs] class PSPAtmosphericalCorrection(ModifierBase): """Correct for atmospheric effects.""" def __call__(self, projectables, optional_datasets=None, **info): """Get the atmospherical correction. Uses pyspectral. """ from pyspectral.atm_correction_ir import AtmosphericalCorrection band = projectables[0] if optional_datasets: satz = optional_datasets[0] else: satz = get_satellite_zenith_angle(band) satz = satz.data # get dask array underneath logger.info("Correction for limb cooling") corrector = AtmosphericalCorrection(band.attrs["platform_name"], band.attrs["sensor"]) atm_corr = da.map_blocks(_call_mapped_correction, satz, band.data, corrector=corrector, band_name=band.attrs["name"], meta=np.array((), dtype=band.dtype)) proj = xr.DataArray(atm_corr, attrs=band.attrs, dims=band.dims, coords=band.coords) self.apply_modifier_info(band, proj) return proj
[docs] class CO2Corrector(ModifierBase): """CO2 correction of the brightness temperature of the MSG 3.9um channel. .. math:: T4_CO2corr = (BT(IR3.9)^4 + Rcorr)^0.25 Rcorr = BT(IR10.8)^4 - (BT(IR10.8)-dt_CO2)^4 dt_CO2 = (BT(IR10.8)-BT(IR13.4))/4.0 Derived from D. Rosenfeld, "CO2 Correction of Brightness Temperature of Channel IR3.9" References: - https://resources.eumetrain.org/IntGuide/PowerPoints/Channels/conversion.ppt """ def __call__(self, projectables, optional_datasets=None, **info): """Apply correction.""" ir_039, ir_108, ir_134 = projectables logger.info("Applying CO2 correction") dt_co2 = (ir_108 - ir_134) / 4.0 rcorr = ir_108 ** 4 - (ir_108 - dt_co2) ** 4 t4_co2corr = (ir_039 ** 4 + rcorr).clip(0.0) ** 0.25 t4_co2corr.attrs = ir_039.attrs.copy() self.apply_modifier_info(ir_039, t4_co2corr) return t4_co2corr