diff --git a/docs/script_examples/plot_weighting_functions.py b/docs/script_examples/plot_weighting_functions.py new file mode 100644 index 00000000..2b919af5 --- /dev/null +++ b/docs/script_examples/plot_weighting_functions.py @@ -0,0 +1,67 @@ +""" +Computation of Weighting Functions +================================== +""" + +# %% +# This example shows how to use the :py:class:`pyrtlib.weighting_functions.WeightingFunctions` method +# to compute the weighting functions for the MWS channels for the U.S. standard atmospheric profile. + +import numpy as np +import warnings +warnings.filterwarnings("ignore", category=UserWarning) +from pyrtlib.weighting_functions import WeightingFunctions +from pyrtlib.climatology import AtmosphericProfiles as atmp +from pyrtlib.utils import ppmv2gkg, mr2rh, get_frequencies_sat + +z, p, _, t, md = atmp.gl_atm(atmp.US_STANDARD) +gkg = ppmv2gkg(md[:, atmp.H2O], atmp.H2O) +rh = mr2rh(p, t, gkg)[0] / 100 + +wf = WeightingFunctions(z, p, t, rh) +wf.frequencies = np.array([50.5, 53.2, 54.35, 54.9, 59.4, 58.825, 58.4]) +wgt = wf.generate_wf() + +wf.plot_wf(wgt, 'Downlooking', ylim=[0, 60], legend=True, figsize=(8, 6), dpi=100) + +#%% +# As above but with the weighting functions computed in uplooking mode. + +wf.satellite = False +wgt = wf.generate_wf() + +wf.plot_wf(wgt, 'Uplooking', ylim=[0, 10], figsize=(8, 6), dpi=100) + +#%% +# The weighting functions can also be computed for a different set of channels. +# The bandpass values are used to compute the weighting functions for the ATMS channels. +# The following code compute the weighting functions for the ATMS channels 5-15. + +cf53 = 53.596 +cf57 = 57.290344 +frq = np.array([52.8, cf53-0.115, cf53+0.115, 54.4, 54.94, 55.5, cf57, + cf57-0.217, cf57+0.217, + cf57-0.3222-0.048, cf57-0.3222+0.048, cf57+0.3222-0.048, cf57+0.3222+0.048, + cf57-0.3222-0.022, cf57-0.3222+0.022, cf57+0.3222-0.022, cf57+0.3222+0.022, + cf57-0.3222-0.010, cf57-0.3222+0.010, cf57+0.3222-0.010, cf57+0.3222+0.010, + cf57-0.3222-0.0045, cf57-0.3222+0.0045, cf57+0.3222-0.0045, cf57+0.3222+0.0045]) + +wf.satellite = True +wf.frequencies = frq +wf.bandpass = np.array([1, 2, 1, 1, 1, 1, 2, 4, 4, 4, 4]) +wf.legend_labels = [f'Channel {i+5}' for i in range(len(wf.bandpass))] +wgt = wf.generate_wf() + +wf.plot_wf(wgt, 'ATMS Channels 5-15', ylim=[0, 70], xlim=[0, 0.11], legend=True, figsize=(8, 6), dpi=100) + +#%% +# The weighting functions can also be computed for a different set of frequencies. +# The following code compute the weighting functions for the MWS channels for a standard tropical atmosphere. +# for grouped frequencies. + +wf.satellite = True +wf.frequencies = get_frequencies_sat('MWS') +wgt = wf.generate_wf() +wf.plot_wf_grouped(wgt, 'MWS Channels (grouped)', ylim=[0, 60], + grouped_frequencies=[4, 9, 19, 1, 13], + grouped_labels=['23-52', '53-55', '57', '89', '164-229'], dpi=350) \ No newline at end of file diff --git a/docs/source/api.rst b/docs/source/api.rst index aa825fb2..f0222245 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -168,6 +168,41 @@ To get all implemented models use the following code: 'R21SD', 'R22SD'] +Weighting Functions +=================== + +Computes the weighting functions to assess the vertical sensitivity of the brightness temperature to the atmospheric profile. + +.. note:: + The weighting functions are computed always using last absorption model available. + +.. autosummary:: + :toctree: generated/ + + pyrtlib.weighting_functions.WeightingFunctions + +.. plot:: + :include-source: true + + from pyrtlib.weighting_functions import WeightingFunctions + from pyrtlib.climatology import AtmosphericProfiles as atmp + from pyrtlib.utils import ppmv2gkg, mr2rh, get_frequencies_sat + + z, p, _, t, md = atmp.gl_atm(atmp.TROPICAL) + gkg = ppmv2gkg(md[:, atmp.H2O], atmp.H2O) + rh = mr2rh(p, t, gkg)[0] / 100 + + wf = WeightingFunctions(z, p, t, rh, .1) + wf.satellite = True + wf.angle = 48. + wf.frequencies = get_frequencies_sat('ICI') + wgt = wf.generate_wf() + + wf.plot_wf_grouped(wgt, '', ylim=[0, 20], + grouped_frequencies=[8, 2, 6, 6, 2], + grouped_labels=['176-190', '240-245', '315-334', '440-455', '659-668']) + + Utility Functions ================= diff --git a/pyrtlib/absorption_model.py b/pyrtlib/absorption_model.py index b8424a30..7aee4053 100644 --- a/pyrtlib/absorption_model.py +++ b/pyrtlib/absorption_model.py @@ -15,6 +15,7 @@ PATH = os.path.dirname(os.path.abspath(__file__)) + class AbsModelError(Exception): """Exception raised for errors in the input model. @@ -119,12 +120,15 @@ def implemented_models() -> Dict[str, List[str]]: 'R22SD'], 'Ozone': ['R18', 'R22']} """ - oxygen_netcdf = Dataset(os.path.join(PATH, '_lineshape', 'o2_lineshape.nc'), mode='r') - wv_netcdf = Dataset(os.path.join(PATH, '_lineshape', 'h2o_lineshape.nc'), mode='r') - ozone_netcdf = Dataset(os.path.join(PATH, '_lineshape', 'o3_lineshape.nc'), mode='r') - - model = {"Oxygen": list(oxygen_netcdf.groups.keys()), - "WaterVapour": list(wv_netcdf.groups.keys()), + oxygen_netcdf = Dataset(os.path.join( + PATH, '_lineshape', 'o2_lineshape.nc'), mode='r') + wv_netcdf = Dataset(os.path.join( + PATH, '_lineshape', 'h2o_lineshape.nc'), mode='r') + ozone_netcdf = Dataset(os.path.join( + PATH, '_lineshape', 'o3_lineshape.nc'), mode='r') + + model = {"Oxygen": list(oxygen_netcdf.groups.keys()), + "WaterVapour": list(wv_netcdf.groups.keys()), "Ozone": list(ozone_netcdf.groups.keys())} return model @@ -279,7 +283,7 @@ def h2o_continuum(self, frq: np.ndarray, vx: np.ndarray, nfreq: int): for j in range(nf): a[j] = 6.532e+12*selfcon[j]*theta**(selftexp[j]+3.) a = np.insert(a, 0, a[1], axis=0) - + for i in range(nfreq): fj = frq/deltaf j = int(fj) @@ -542,12 +546,12 @@ def h2o_absorption(self, pdrykpa: np.ndarray, vx: np.ndarray, ekpa: np.ndarray, # separate the following original equ. into line and continuum # terms, and change the units from np/km to ppm # abh2o = .3183e-4*den*sum + con - + if H2OAbsModel.model in ['R23SD', 'R24']: conf = self.h2oll.cf * ti**self.h2oll.xcf con = (conf * pda + npp_cs * pvap) * pvap * f**2 - h20m = 2.9915075E-23 # mass of water molecule (g) + h20m = 2.9915075E-23 # mass of water molecule (g) if H2OAbsModel.model in ['R22SD', 'R23SD']: npp = 1.e-10 * rho * summ / (np.pi * h20m) / db2np / factor elif H2OAbsModel.model == 'R24': @@ -699,7 +703,7 @@ def o2_absorption(self, pdrykpa: float, vx: float, ekpa: float, frq: float, amu: pe1 = .99 * den pe2 = pe1**2 if O2AbsModel.model in ['R24']: - pe1 = den # [8] includes common TEMP-dependence + pe1 = den # [8] includes common TEMP-dependence pe2 = pe1**2 dfnr = self.o2ll.wb300 * den @@ -741,13 +745,14 @@ def o2_absorption(self, pdrykpa: float, vx: float, ekpa: float, frq: float, amu: a = np.zeros(nlines) g = np.zeros(nlines) for k in range(0, nlines): - a[k] = self.o2ll.s300[k]*np.exp(-self.o2ll.be[k] * th1)/self.o2ll.f[k]**2 + a[k] = self.o2ll.s300[k] * \ + np.exp(-self.o2ll.be[k] * th1)/self.o2ll.f[k]**2 g[k] = self.o2ll.g0[k] + self.o2ll.g1[k]*th1 if k > 0 and k <= 37: summ += a[k]*g[k] anorm += a[k] off = summ/anorm - summ = (1.584e-17/th) * dfnr/ (freq * freq + dfnr * dfnr) + summ = (1.584e-17/th) * dfnr / (freq * freq + dfnr * dfnr) for k in range(0, nlines): width = self.o2ll.w300[k] * den y = pe1*(self.o2ll.y300[k]+self.o2ll.y1[k]*th1) @@ -755,17 +760,19 @@ def o2_absorption(self, pdrykpa: float, vx: float, ekpa: float, frq: float, amu: gfac = 1. + pe2 * (g[k] - off) else: gfac = 1. - fcen = self.o2ll.f[k] + pe2 * (self.o2ll.dnu0[k] + self.o2ll.dnu1[k] * th1) + fcen = self.o2ll.f[k] + pe2 * \ + (self.o2ll.dnu0[k] + self.o2ll.dnu1[k] * th1) if k == 0 and np.abs(freq-fcen) < 10. * width: width2 = .076 * width xc = complex(width-1.5*width2, freq-fcen)/width2 xrt = np.sqrt(xc) pxw = 1.77245385090551603 * xrt * \ - _dcerror(-np.imag(xrt), np.real(xrt)) - asd = complex(1.,y)*2.*(1.-pxw)/width2 + _dcerror(-np.imag(xrt), np.real(xrt)) + asd = complex(1., y)*2.*(1.-pxw)/width2 sf1 = np.real(asd) else: - sf1 = (width*gfac + (freq-fcen)*y)/((freq-fcen)**2 + width**2) + sf1 = (width*gfac + (freq-fcen)*y) / \ + ((freq-fcen)**2 + width**2) sf2 = (width*gfac - (freq+fcen)*y)/((freq+fcen)**2 + width**2) summ += a[k] * (sf1+sf2) if k == 37: @@ -776,15 +783,17 @@ def o2_absorption(self, pdrykpa: float, vx: float, ekpa: float, frq: float, amu: sumy = 1.584e-17 * self.o2ll.wb300 sumg = 0. asq = 0. - adjy = .99 # adjustment following update of line intensities + adjy = .99 # adjustment following update of line intensities a = np.zeros(nlines) y = np.zeros(nlines) g = np.zeros(nlines) for k in range(0, nlines): - a[k] = self.o2ll.s300[k] * np.exp(-self.o2ll.be[k] * th1) * th/self.o2ll.f[k]**2 + a[k] = self.o2ll.s300[k] * \ + np.exp(-self.o2ll.be[k] * th1) * th/self.o2ll.f[k]**2 y[k] = adjy * (self.o2ll.y300[k] + self.o2ll.y1[k]*th1) if k <= 37: - sumy += 2. * a[k] * (self.o2ll.w300[k] + y[k] * self.o2ll.f[k]) + sumy += 2. * a[k] * \ + (self.o2ll.w300[k] + y[k] * self.o2ll.f[k]) g[k] = self.o2ll.g0[k] + self.o2ll.g1[k] * th1 if k > 0 and k <= 37: sumg += a[k] * g[k] @@ -795,9 +804,9 @@ def o2_absorption(self, pdrykpa: float, vx: float, ekpa: float, frq: float, amu: sumy2 = sumy/(2. * anorm) ratio = sumg/asq for k in range(1, 38): - y[k] -= sumy2/self.o2ll.f[k] # bias adjustment - g[k] -= a[k] * ratio # makes G orthogonal to A - + y[k] -= sumy2/self.o2ll.f[k] # bias adjustment + g[k] -= a[k] * ratio # makes G orthogonal to A + summ = 1.584e-17 * wnr/(freq * freq + wnr * wnr) for k in range(0, nlines): width = self.o2ll.w300[k] * pe1 @@ -806,17 +815,19 @@ def o2_absorption(self, pdrykpa: float, vx: float, ekpa: float, frq: float, amu: gfac = 1. + pe2 * g[k] else: gfac = 1. - fcen = self.o2ll.f[k] + pe2 * (self.o2ll.dnu0[k] + self.o2ll.dnu1[k] * th1) + fcen = self.o2ll.f[k] + pe2 * \ + (self.o2ll.dnu0[k] + self.o2ll.dnu1[k] * th1) if k == 0 and np.abs(freq-fcen) < 10. * width: width2 = .076 * width xc = complex(width - 1.5*width2, (freq-fcen))/width2 xrt = np.sqrt(xc) pxw = 1.77245385090551603 * xrt * \ - _dcerror(-np.imag(xrt), np.real(xrt)) - asd = complex(1.,yk) * 2. * (1.-pxw)/width2 + _dcerror(-np.imag(xrt), np.real(xrt)) + asd = complex(1., yk) * 2. * (1.-pxw)/width2 sf1 = np.real(asd) else: - sf1 = (width*gfac + (freq-fcen)*yk)/((freq-fcen)**2 + width**2) + sf1 = (width*gfac + (freq-fcen)*yk) / \ + ((freq-fcen)**2 + width**2) sf2 = (width*gfac - (freq+fcen)*yk)/((freq+fcen)**2 + width**2) summ += a[k] * (sf1+sf2) else: diff --git a/pyrtlib/rt_equation.py b/pyrtlib/rt_equation.py index d6341448..d18425f9 100644 --- a/pyrtlib/rt_equation.py +++ b/pyrtlib/rt_equation.py @@ -357,7 +357,8 @@ def cloud_integrated_density(dencld: np.ndarray, ds: np.ndarray, lbase: np.ndarr scld = 0.0 for i in range(0, ncld): for j in range(int(lbase[i]) + 1, int(ltop[i])): - scld += np.dot(ds[j], (np.dot(0.5, (dencld[j] + dencld[j - 1])))) + scld += np.dot(ds[j], + (np.dot(0.5, (dencld[j] + dencld[j - 1])))) # convert the integrated value to cm. scld = np.dot(scld, 0.1) @@ -562,7 +563,8 @@ def clearsky_absorption(p: np.ndarray, t: np.ndarray, e: np.ndarray, frq: np.nda 'No model avalaible with this name: {} . Sorry...'.format('model')) if isinstance(o3n, np.ndarray) and O3AbsModel.model in ['R18', 'R21', 'R21SD', 'R22', 'R22SD', 'R23', 'R23SD', 'R24']: - aO3[i] = O3AbsModel().o3_absorption(t[i], p[i], frq, o3n[i], amu) + aO3[i] = O3AbsModel().o3_absorption( + t[i], p[i], frq, o3n[i], amu) adry[i] = aO2[i] + aN2[i] + aO3[i] diff --git a/pyrtlib/tb_spectrum.py b/pyrtlib/tb_spectrum.py index 088270b7..5ec702cc 100644 --- a/pyrtlib/tb_spectrum.py +++ b/pyrtlib/tb_spectrum.py @@ -106,6 +106,8 @@ def __init__(self, z: np.ndarray, p: np.ndarray, t: np.ndarray, rh: np.ndarray, self.sptauwet = np.zeros((self.nf, self.nang)) self.sptauliq = np.zeros((self.nf, self.nang)) self.sptauice = np.zeros((self.nf, self.nang)) + self.awet = np.zeros((self.nf, self.nang, self.nl)) + self.adry = np.zeros((self.nf, self.nang, self.nl)) self.ptaudry = np.zeros((self.nf, self.nang, self.nl)) self.ptaulay = np.zeros((self.nf, self.nang, self.nl)) self.ptauwet = np.zeros((self.nf, self.nang, self.nl)) @@ -345,14 +347,14 @@ def execute(self, only_bt: bool = True) -> Union[pd.DataFrame, Tuple[pd.DataFram for j in range(0, self.nf): RTEquation._emissivity = self._es[j] # Rosenkranz, personal communication, 2019/02/12 (email) - awet, adry = RTEquation.clearsky_absorption(self.p, self.tk, e, self.frq[j], + self.awet[j, k, :], self.adry[j, k, :] = RTEquation.clearsky_absorption(self.p, self.tk, e, self.frq[j], self.o3n, self.amu if self._uncertainty else None) self.sptauwet[j, k], \ self.ptauwet[j, k, :] = RTEquation.exponential_integration( - True, awet, ds, 1, self.nl, 1) + True, self.awet[j, k, :], ds, 1, self.nl, 1) self.sptaudry[j, k], \ self.ptaudry[j, k, :] = RTEquation.exponential_integration( - True, adry, ds, 1, self.nl, 1) + True, self.adry[j, k, :], ds, 1, self.nl, 1) if self.cloudy: aliq, aice = RTEquation.cloudy_absorption( self.tk, self.denliq, self.denice, self.frq[j]) @@ -398,5 +400,6 @@ def execute(self, only_bt: bool = True) -> Union[pd.DataFrame, Tuple[pd.DataFram else: return df, {'taulaywet': self.ptauwet, 'taulaydry': self.ptaudry, 'taulayliq': self.ptauliq, 'taulayice': self.ptauice, + 'awet': self.awet, 'adry': self.adry, 'srho': self.srho, 'swet': self.swet, 'sdry': self.sdry, 'sliq': self.sliq, 'sice': self.sice} diff --git a/pyrtlib/tests/test_wf.py b/pyrtlib/tests/test_wf.py new file mode 100644 index 00000000..a45f76af --- /dev/null +++ b/pyrtlib/tests/test_wf.py @@ -0,0 +1,57 @@ +from unittest import TestCase + +import numpy as np +from pyrtlib.weighting_functions import WeightingFunctions +from pyrtlib.climatology import AtmosphericProfiles as atmp +from pyrtlib.utils import ppmv2gkg, mr2rh +import numpy as np +import matplotlib +matplotlib.use('Template') + +class Test(TestCase): + def wf_computation(self, frq: np.ndarray): + z, p, d, t, md = atmp.gl_atm(atmp.TROPICAL) + gkg = ppmv2gkg(md[:, atmp.H2O], atmp.H2O) + rh = mr2rh(p, t, gkg)[0] / 100 + + wf = WeightingFunctions(z, p, t, rh) + wf.frequencies = frq + + return wf.generate_wf(), z, wf + + def test_f_89(self): + wgt, z, _ = self.wf_computation(np.array([89.0])) + np.testing.assert_almost_equal(z[np.argmax(wgt)], 0., decimal=15) + + def test_f_57(self): + wgt, z, _ = self.wf_computation(np.array([57.290344])) + np.testing.assert_almost_equal(z[np.argmax(wgt)], 18., decimal=15) + + def test_f_57_6(self): + wgt, z, _ = self.wf_computation(np.array([57.660544])) + np.testing.assert_almost_equal(z[np.argmax(wgt)], 27.5, decimal=15) + + def test_plot_wf(self): + wgt, z, wf = self.wf_computation(np.array([57.660544])) + wf.plot_wf(wgt, 'Test', legend=True, ylim=(0, 60), xlim=(0, .1), figsize=(8, 6), dpi=100) + + def test_plot_wf_grouped(self): + wgt, z, wf = self.wf_computation(np.array([57.660544])) + wf.plot_wf_grouped(wgt, 'Test', grouped_frequencies=[ + 1], grouped_labels=['Test'], dpi=100) + + def test_plot_wf_multiple(self): + freqs = np.array([57.290344, 57.660544, 89.0]) + wgt, z, wf = self.wf_computation(freqs) + wf.satellite = False + wf.plot_wf(wgt, 'Multiple Frequencies', + legend=True, figsize=(8, 6), dpi=100) + + def test_plot_wf_bandpass(self): + cf53 = 53.596 + cf57 = 57.290344 + freqs = np.array([52.8, cf53-0.115, cf53+0.115, 54.4, 54.94, 55.5, cf57, cf57-0.217, cf57+0.217, cf57-0.3222-0.048, cf57-0.3222+0.048, cf57+0.3222-0.048, cf57+0.3222+0.048]) + wgt, z, wf = self.wf_computation(freqs) + wf.bandpass = np.array([1, 2, 1, 1, 1, 1, 2, 4]) + wf.plot_wf(wgt, 'Bandpass', + legend=True, figsize=(8, 6), dpi=100) diff --git a/pyrtlib/utils.py b/pyrtlib/utils.py index 66892b4c..2e6682ad 100644 --- a/pyrtlib/utils.py +++ b/pyrtlib/utils.py @@ -323,7 +323,8 @@ def mr2e(p: np.ndarray, mr: np.ndarray) -> np.ndarray: return e -def rho2rh(rho: np.ndarray, t: np.ndarray, p: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: + +def rho2rh(rho: np.ndarray, t: np.ndarray, p: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: """Convert water vapor density to relative humidity. Args: @@ -908,11 +909,12 @@ def get_frequencies(instrument: Optional[str] = 'hat') -> List: ], 'k2w': [23.8400, 31.4000, 72.5000, 82.5000, 90.0000, 150.000] } - + try: return np.array(frequencies[instrument.lower()]) except KeyError: - raise ValueError(f"Invalid instrument name. Available instruments are: {list(frequencies.keys())}") + raise ValueError( + f"Invalid instrument name. Available instruments are: {list(frequencies.keys())}") def to_kelvin(t: np.ndarray) -> np.ndarray: @@ -970,7 +972,7 @@ def get_frequencies_sat(instrument: str) -> np.ndarray: cf165 = 165.5 cf53 = 57.290344 cf57 = 57.290344 - + instrument = instrument.upper() if instrument == 'SAPHIR': diff --git a/pyrtlib/weighting_functions.py b/pyrtlib/weighting_functions.py new file mode 100644 index 00000000..ba709e07 --- /dev/null +++ b/pyrtlib/weighting_functions.py @@ -0,0 +1,377 @@ +# -*- coding: utf-8 -*- +""" +This module is used to define the weighting functions +""" + +__author__ = '' +__date__ = 'May 2024' +__copyright__ = '(C) 2024, CNR-IMAA' + +from typing import Dict, Tuple, Optional, List +import warnings +# warnings.filterwarnings("ignore", category=UserWarning) +import numpy as np +import matplotlib.pyplot as plt +from . absorption_model import AbsModel +from .tb_spectrum import TbCloudRTE +from .absorption_model import O2AbsModel + + +class WeightingFunctions(object): + """ + This class is used to compute the weighting functions + + .. note:: + The weighting functions are computed always using last absorption model available. + """ + + def __init__(self, z: np.ndarray, p: np.ndarray, t: np.ndarray, rh: np.ndarray, + value_prof_interpolation: Optional[float] = 1) -> None: + """ + Constructor for the WeightingFunctions class. + + Args: + z (np.ndarray): Height profile. + p (np.ndarray): Pressure profile. + t (np.ndarray): Temperature profile. + rh (np.ndarray): Relative humidity profile. + value_prof_interpolation (Optional[float], optional): The value for the interpolation of the profiles. Defaults to 1. + + Raises: + ValueError: If the input values are not valid arrays. + + Returns: + None + + .. note:: + The input arrays must have the same size. Units for the arrays are: + - z: km + - p: hPa + - t: K + - rh: fraction + - value_prof_interpolation: km + """ + + self.zinterp = np.arange( + 0, np.max(z)+value_prof_interpolation, value_prof_interpolation) + self.pinterp = np.interp(self.zinterp, z, p) + self.tinterp = np.interp(self.zinterp, z, t) + self.rhinterp = np.interp(self.zinterp, z, rh) + + self._frequencies = None + self._satellite = True + self.emissivity = 1. + self.model = AbsModel.implemented_models()['WaterVapour'][-1] + self.angle = 90. + self.bandpass = None + self.legend_labels = None + + @property + def satellite(self) -> bool: + """If :code:`True` computes an upward-propagating brightness-temperature spectrum + otherwise a downward-propagating brightness-temperature + spectrum at the bottom of the atmosphere will be performed. + """ + return self._satellite + + @satellite.setter + def satellite(self, sat: bool) -> None: + """ + Sets the satellite flag. + + Parameters: + - sat (bool): A boolean value indicating whether the satellite flag should be set to True or False. + + Raises: + - ValueError: If the input value for `sat` is not a boolean. + + Returns: + - None + """ + if isinstance(sat, bool): + self._satellite = sat + else: + raise ValueError("Please enter True or False") + + @property + def frequencies(self) -> np.ndarray: + """Returns the frequencies array. + + Returns: + np.ndarray: The frequencies array. + """ + return self._frequencies + + @frequencies.setter + def frequencies(self, frequencies: np.ndarray) -> None: + """ + Set the frequencies for the weighting function. + + Parameters: + - frequencies: np.ndarray + The frequencies to be set. + + Raises: + - ValueError: If the input is not a valid array for frequencies. + + Returns: + - None + """ + if isinstance(frequencies, np.ndarray): + self._frequencies = frequencies + else: + raise ValueError( + "Please enter a valid array for frequencies") + + def _perform_calculation(self) -> Dict[str, np.ndarray]: + """ + Perform the calculation to get opacity and absorption. + """ + + rte = TbCloudRTE(self.zinterp, self.pinterp, self.tinterp, + self.rhinterp, self.frequencies, angles=np.array([self.angle]), ray_tracing=True) + rte.satellite = self.satellite + rte.init_absmdl(self.model) + rte.emissivity = self.emissivity + O2AbsModel.model = AbsModel.implemented_models()['Oxygen'][-1] + O2AbsModel.set_ll() + _, other = rte.execute(False) + + return other + + def _mean_channel(self, a: np.ndarray, bandpass: np.ndarray) -> np.ndarray: + """ + Compute the mean of the input array based on bandpass. + + Args: + a (np.ndarray): The input array. + bandpass (np.ndarray): The bandpass array. + + Returns: + np.ndarray: The mean of the input array based on bandpass. + """ + + channel_mean = np.zeros((len(bandpass), a.shape[1])) + cnt = 0 + for i, _ in enumerate(bandpass): + channel_mean[i, :] = np.mean(a[cnt:bandpass[i]+cnt, :], axis=0) + cnt += bandpass[i] + + return channel_mean + + def generate_wf(self) -> np.ndarray: + r""" + Generate the weighting function. + + .. math:: \frac{\partial \tau}{\partial s} = -{\tau (\nu, s)}\times{\alpha (\nu, s)} + + This method calculates the weighting function based on the provided data. + It performs calculations and returns the resulting weighting function. + + Returns: + np.ndarray: The calculated weighting function. + + """ + other = self._perform_calculation() + if self.satellite: + z0 = self.zinterp[-1] + else: + z0 = self.zinterp[0] + self.l0 = np.argmin(np.abs(self.zinterp-z0)) + + tautotal = (other['taulaywet'] + other['taulaydry'] + + other['taulayliq'] + other['taulayice']) + abstotal = other['awet'] + other['adry'] + + if self.satellite: + ao = np.fliplr(tautotal[:, 0, :self.l0+1]) + abss = abstotal[:, 0, :self.l0+1] + else: + ao = tautotal[:, 0, self.l0:] + abss = abstotal[:, 0, self.l0:] + + if self.satellite: + cumao = np.fliplr(np.cumsum(ao[::-1], axis=1)[::-1]) + else: + cumao = np.cumsum(ao, axis=1) + trans = np.exp(-cumao) + wgt = np.array([trans[:, j] * abss[:, j] + for j in range(trans.shape[1])]).T + + return wgt + + def plot_wf(self, wgt: np.ndarray, title: Optional[str] = '', ylim: Optional[Tuple[int, int]] = None, + xlim: Optional[Tuple[int, int]] = None, + legend: Optional[bool] = True, + normalized: Optional[bool] = False, + **kwargs) -> None: + """ + Plot the weighting functions + + Args: + wgt (ndarray): The weighting functions to be plotted. + title (Optional[str], optional): The title of the plot. Defaults to ''. + ylim (Optional[Tuple[int, int]], optional): The y-axis limits of the plot. Defaults to None. + xlim (Optional[Tuple[int, int]], optional): The x-axis limits of the plot. Defaults to None. + legend (Optional[bool], optional): Whether to show the legend. Defaults to True. + normalized (Optional[bool], optional): Whether to normalize to the max the weighting functions. Defaults to False. + **kwargs: Additional keyword arguments (figsize, dpi) + + """ + plt.rcdefaults() + plt.rcParams.update({'font.size': 15}) + plt.rcParams['font.family'] = 'Arial' + plt.rcParams['font.stretch'] = 'condensed' + # plt.rcParams["font.weight"] = "bold" + # plt.rcParams["axes.labelweight"] = "bold" + # plt.rc('font', family=['cmr10']) + # plt.rc('mathtext', fontset='cm') + # plt.rc('axes.formatter', use_mathtext=False) + + if self.bandpass is not None: + wgt = self._mean_channel(wgt, self.bandpass) + self.frequencies = np.array( + [np.mean(self.frequencies[i:i+j]) for i, j in enumerate(self.bandpass)]) + + # if yaxis == 'p': + # y = self.pinterp[:-1] + # y_label = 'Pressure [$hPa$]' + # units = 'hPa' + # else: + y = self.zinterp[:-1] + y_label = 'Height [$km$]' + units = 'km' + + normalize = np.max( + wgt, axis=1) if normalized else np.ones(wgt.shape[1]) + + plt.figure(figsize=kwargs.get('figsize', (5, 10)), + dpi=kwargs.get('dpi', 150)) + for i, frequenciy in enumerate(self.frequencies): + if self.satellite: + line_s = 'dashed' if y[np.argmax(wgt[i, 1:])] == 0. else None + plt.plot(wgt[i, 1:]/normalize[i], y, linestyle=line_s, + label=f'{frequenciy:.2f} GHz {y[np.argmax(wgt[i, 1:])]:.0f} {units}') + else: + plt.plot(wgt[i, :], self.zinterp, + label=f'{frequenciy:.2f} GHz') + + # if yaxis == 'p': + # plt.gca().invert_yaxis() + # plt.yscale('log') + plt.xlabel('Weighting Function [$km^{-1}$]') + plt.ylabel(y_label) + if legend: + if self.legend_labels: + plt.legend(self.legend_labels, fontsize=10) + else: + plt.legend(fontsize=10) + if ylim: + plt.ylim(ylim[0], ylim[1]) + if xlim: + plt.xlim(xlim[0], xlim[1]) + + if title: + plt.title(title) + plt.grid(color='#d6d6d6', linestyle='dashed') + plt.tick_params(axis='both', which='major', direction='in') + ax_twinx = plt.twinx() + ax_twinx.set_ylim((self.pinterp[np.where(self.zinterp == np.min(self.zinterp))], + self.pinterp[np.where(self.zinterp == np.max(self.zinterp))])) + if ylim: + ax_twinx.set_ylim((self.pinterp[np.where(self.zinterp == ylim[0])], + self.pinterp[np.where(self.zinterp == ylim[1])])) + ax_twinx.set_yscale('log') + ax_twinx.set_ylabel('Pressure [$hPa$]') + ax_twinx.tick_params(axis='both', which='both', direction='in') + plt.tight_layout() + plt.show() + + def plot_wf_grouped(self, wgt: np.ndarray, title: Optional[str] = '', + ylim: Optional[Tuple[int, int]] = None, + grouped_frequencies: List[int] = None, + grouped_labels: List[str] = None, + **kwargs) -> None: + """ + Plot the weighting functions + + Args: + wgt (ndarray): The weighting functions to be plotted. + title (Optional[str], optional): The title of the plot. Defaults to ''. + ylim (Optional[Tuple[int, int]], optional): The y-axis limits of the plot. Defaults to None. + grouped_frequencies (List[int]): The number of frequencies to be grouped. + grouped_labels (List[str]): The labels for the grouped frequencies. + dpi (Optional[int], optional): The dpi of the plot. Defaults to 150. + **kwargs: Additional keyword arguments (figsize, dpi) + + Raises: + ValueError: If the plot is not available for the current satellite mode. + """ + if not self.satellite: + raise ValueError("This plot is only available for satellite mode") + + plt.rcdefaults() + plt.rcParams.update({'font.size': 15}) + # plt.rc('font', family=['Sans Serif']) + plt.rcParams['font.family'] = 'Arial' + plt.rcParams['font.stretch'] = 'condensed' + # plt.rcParams["font.weight"] = "bold" + # plt.rcParams["axes.labelweight"] = "bold" + # plt.rc('mathtext', fontset='cm') + # plt.rc('axes.formatter', use_mathtext=False) + + y = self.zinterp[:-1] + y_label = 'Height [$km$]' + + plt.figure(figsize=kwargs.get('figsize', (10, 4)), + dpi=kwargs.get('dpi', 150)) + cnt = 0 + max_increment = np.array([]) + for i in grouped_frequencies: + a = np.tile(cnt, i) + max_increment = np.append(max_increment, a) + cnt += 0.40 + + cnt = 0 + g_labels = np.array([], dtype='object') + for i, grouped_frequency in enumerate(grouped_frequencies): + a = np.repeat(grouped_labels[i], grouped_frequency) + g_labels = np.append(g_labels, a) + cnt += 1 + + for i in range(len(self.frequencies)): + line_s = 'dashed' if y[np.argmax(wgt[i, 1:])] == 0. else None + plt.plot(wgt[i, 1:]+max_increment[i], y, + linestyle=line_s, linewidth=1.5) + + for i in range(len(grouped_frequencies)): + index_wgt = np.cumsum(np.array(grouped_frequencies)) - 1 + height = ylim[1]-1 if ylim else np.max(y)-1 + perc_height = height/15 + incr_wgt = np.min(wgt[index_wgt[i], 1:]) + \ + max_increment[index_wgt[i]] + plt.text(.02 + incr_wgt, height-perc_height, + f'{grouped_labels[i]}', ha='left', va='top') + + plt.xlabel('Weighting Function [$km^{-1}$]') + plt.ylabel(y_label) + if ylim: + plt.ylim(ylim[0], ylim[1]) + if title: + plt.title(title) + plt.grid(color='#d6d6d6') + plt.tick_params(axis='both', which='major', direction='in') + ax_twinx = plt.twinx() + ax_twinx.set_ylim((self.pinterp[np.where(self.zinterp == np.min(self.zinterp))], + self.pinterp[np.where(self.zinterp == np.max(self.zinterp))])) + if ylim: + ax_twinx.set_ylim((self.pinterp[np.where(self.zinterp == ylim[0])], + self.pinterp[np.where(self.zinterp == ylim[1])])) + ax_twinx.set_yscale('log') + ax_twinx.set_ylabel('Pressure [$hPa$]') + ax_twinx.tick_params(axis='both', which='both', direction='in') + plt.tight_layout() + plt.show() + + warnings.warn( + "Weighting functions have been shifted to the right so that they all can be displayed in a single graph") diff --git a/requirements.txt b/requirements.txt index bd62b0e8..3cf6d320 100644 --- a/requirements.txt +++ b/requirements.txt @@ -6,3 +6,4 @@ scikit-learn netCDF4 requests bs4 +matplotlib