Merge pull request #26 from SatCloP/wgt

[feature] add computation of weighting functions

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)

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
 =================
 

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,31 +745,34 @@ 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)
                 if k > 0 and k <= 37:
                     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:

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]
 

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}

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)