add subclass in Result for Active and Passive mode to avoid 'beginner…

…' mistake to call Tb for an Active sensor and sigma for a Passive sensor.

smrt/core/result.py

@@ -1,16 +1,19 @@
 # coding: utf-8
 
 """ The results of RT Solver are hold by the :py:class:`Result` class. This class provides several functions
-to access to the Stokes Vector and Muller matrix in a simple way. Most notable ones are :py:meth:`Result.TbV` and :py:meth:`Result.TbH` for the passive mode calculations and
-:py:meth:`Result.sigmaHH` and :py:meth:`Result.sigmaVV`. Other methods could be developed for specific needs.
+to access to the Stokes Vector and Muller matrix in a simple way. Most notable ones are :py:meth:`Result.TbV` and :py:meth:`Result.TbH`
+for the passive mode calculations and :py:meth:`Result.sigmaHH` and :py:meth:`Result.sigmaVV`. Other methods could be developed for
+specific needs.
 
-To save results of calculations in a file, simply use the pickle module or other serialization schemes. We may provide a unified and inter-operable solution in the future.
+To save results of calculations in a file, simply use the pickle module or other serialization schemes. We may provide a unified and
+inter-operable solution in the future.
 
-Under the hood, :py:class:`Result` uses xarray module which provides multi-dimensional array with explicit, named, dimensions. Here the common dimensions
-are frequency, polarization, polarization_inc, theta_inc, theta, and phi. They are created by the RT Solver. The interest of using named dimension is that slice of the xarray (i.e. results)
-can be selected based on the dimension name whereas with numpy the order of the dimensions matters. Because this is very convenient, users may be
-interested in adding other dimensions specific to their context such as time, longitude, latitude, points, ... To do so, :py:meth:`smrt.core.model.Model.run` accepts a list of snowpack
-and optionally the parameter `snowpack_dimension` is used to specify the name and values of the new dimension to build.
+Under the hood, :py:class:`Result` uses xarray module which provides multi-dimensional array with explicit, named, dimensions. Here the
+common dimensions are frequency, polarization, polarization_inc, theta_inc, theta, and phi. They are created by the RT Solver. The interest
+ of using named dimension is that slice of the xarray (i.e. results) can be selected based on the dimension name whereas with numpy the order
+ of the dimensions matters. Because this is very convenient, users may be interested in adding other dimensions specific to their context such
+  as time, longitude, latitude, points, ... To do so, :py:meth:`smrt.core.model.Model.run` accepts a list of snowpack and optionally
+  the parameter `snowpack_dimension` is used to specify the name and values of the new dimension to build.
 
 Example::
 
@@ -33,6 +36,7 @@
 import xarray as xr
 import pandas as pd
 from smrt.utils import dB
+from smrt.core.error import SMRTError
 
 
 def open_result(filename):
@@ -44,7 +48,23 @@ def open_result(filename):
         if d.startswith("polarization"):
             data[d] = data[d].astype("U1")
 
-    return Result(data)
+    mode = getattr(data.attrs, 'mode', None)
+    if (mode is None) or (mode not in 'AP'):
+        # guess the mode
+        if 'theta_inc' in data.coords:
+            mode = 'A'
+        else:
+            mode = 'P'
+
+    return make_result(mode, data)
+
+
+def make_result(mode, *args, **kwargs):
+    """create an active or passive result object according to the mode"""
+    if mode == 'A':
+        return ActiveResult(*args, **kwargs)
+    else:
+        return PassiveResult(*args, **kwargs)
 
 
 class Result(object):
@@ -61,96 +81,130 @@ def __init__(self, intensity, coords=None):
         else:
             self.data = xr.DataArray(intensity, coords)
 
+        if hasattr(self, "mode"):
+            self.data.attrs['mode'] = self.mode
+        else:
+            raise SMRTError("Result base class is abstract, uses a subclass instead. The subclass must define the 'mode' attribute")
+
     @property
     def coords(self):
-        """Return the coordinates of the result (theta, frequency, ...). Note that the coordinates are also result attribute, so result.frequency works (and so on for all the coordinates)."""
+        """Return the coordinates of the result (theta, frequency, ...). Note that the coordinates are also result attribute,
+        so result.frequency works (and so on for all the coordinates)."""
         return self.data.coords
 
     def __getattr__(self, attr):
         if attr != "data" and attr in self.data.coords:
             return self.data.coords[attr]
         else:
-            return super().__getattr__(attr)
+            raise AttributeError("AttributeError: '%s' object has no attribute '%s'" % (type(self), attr))
+
+    def save(self, filename):
+        """save a result to disk. Under the hood, this is a netCDF file produced by xarray (http://xarray.pydata.org/en/stable/io.html)."""
+        self.data.to_netcdf(filename)
+
+
+class PassiveResult(Result):
+
+    mode = 'P'
 
     def Tb(self, **kwargs):
-        """Return brightness temperature. Any parameter can be added to slice the results (e.g. frequency=37e9 or polarization='V'). See xarray slicing with sel method (to document)"""
+        """Return brightness temperature. Any parameter can be added to slice the results (e.g. frequency=37e9 or polarization='V').
+         See xarray slicing with sel method (to document)"""
         return _strongsqueeze(self.data.sel(**kwargs))
 
     def Tb_as_dataframe(self, **kwargs):
-        """Return brightness temperature. Any parameter can be added to slice the results (e.g. frequency=37e9 or polarization='V'). See xarray slicing with sel method (to document)"""
+        """Return brightness temperature. Any parameter can be added to slice the results (e.g. frequency=37e9 or polarization='V').
+         See xarray slicing with sel method (to document)"""
         return self.Tb(**kwargs).to_dataframe(name='Tb')
 
     def TbV(self, **kwargs):
-        """Return V polarization. Any parameter can be added to slice the results (e.g. frequency=37e9). See xarray slicing with sel method (to document)"""
+        """Return V polarization. Any parameter can be added to slice the results (e.g. frequency=37e9).
+         See xarray slicing with sel method (to document)"""
         return _strongsqueeze(self.data.sel(polarization='V', **kwargs))
 
     def TbH(self, **kwargs):
-        """Return H polarization. Any parameter can be added to slice the results (e.g. frequency=37e9). See xarray slicing with sel method (to document)"""
+        """Return H polarization. Any parameter can be added to slice the results (e.g. frequency=37e9).
+         See xarray slicing with sel method (to document)"""
         return _strongsqueeze(self.data.sel(polarization='H', **kwargs))
 
     def polarization_ratio(self, ratio="H_V", **kwargs):
-        """Return polarization ratio. Any parameter can be added to slice the results (e.g. frequency=37e9). See xarray slicing with sel method (to document)"""
+        """Return polarization ratio. Any parameter can be added to slice the results (e.g. frequency=37e9).
+         See xarray slicing with sel method (to document)"""
         return _strongsqueeze(self.data.sel(polarization=ratio[0], **kwargs) / self.data.sel(polarization=ratio[-1], **kwargs))
 
+
+class ActiveResult(Result):
+
+    mode = 'A'
+
     def sigma(self, **kwargs):
-        """Return backscattering coefficient. Any parameter can be added to slice the results (e.g. frequency=37e9, polarization_inc='V', polarization='V'). See xarray slicing with sel method (to document)"""
+        """Return backscattering coefficient. Any parameter can be added to slice the results (e.g. frequency=37e9, polarization_inc='V', polarization='V').
+         See xarray slicing with sel method (to document)"""
         if 'theta' in kwargs:
             theta = np.atleast_1d(kwargs.pop('theta'))
         else:
             theta = self.data.theta
 
         backscatter = self.data.sel(**kwargs).sel_points(theta_inc=theta, theta=theta).rename({'points': 'theta'})
-        return 4*np.pi*_strongsqueeze(np.cos(np.radians(theta)) * backscatter)
+        return 4 * np.pi * _strongsqueeze(np.cos(np.radians(theta)) * backscatter)
 
     def sigma_dB(self, **kwargs):
-        """Return backscattering coefficient. Any parameter can be added to slice the results (e.g. frequency=37e9, polarization_inc='V', polarization='V'). See xarray slicing with sel method (to document)"""
+        """Return backscattering coefficient. Any parameter can be added to slice the results (e.g. frequency=37e9, polarization_inc='V', polarization='V').
+         See xarray slicing with sel method (to document)"""
         return dB(self.simga(**kwargs))
 
     def sigma_as_dataframe(self, **kwargs):
-        """Return backscattering coefficient. Any parameter can be added to slice the results (e.g. frequency=37e9, polarization_inc='V', polarization='V'). See xarray slicing with sel method (to document)"""
+        """Return backscattering coefficient. Any parameter can be added to slice the results (e.g. frequency=37e9, polarization_inc='V', polarization='V').
+         See xarray slicing with sel method (to document)"""
         return self.sigma(**kwargs).to_dataframe(name='sigma')
 
     def sigma_dB_as_dataframe(self, **kwargs):
-        """Return backscattering coefficient in dB. Any parameter can be added to slice the results (e.g. frequency=37e9, polarization_inc='V', polarization='V'). See xarray slicing with sel method (to document)"""
+        """Return backscattering coefficient in dB. Any parameter can be added to slice the results (e.g. frequency=37e9, polarization_inc='V', 
+        polarization='V'). See xarray slicing with sel method (to document)"""
         return self.sigma_dB(**kwargs).to_dataframe(name='sigma')
 
     def sigmaVV(self, **kwargs):
-        """Return VV backscattering coefficient. Any parameter can be added to slice the results (e.g. frequency=37e9). See xarray slicing with sel method (to document)"""
+        """Return VV backscattering coefficient. Any parameter can be added to slice the results (e.g. frequency=37e9).
+         See xarray slicing with sel method (to document)"""
         return self.sigma(polarization_inc='V', polarization='V', **kwargs)
 
     def sigmaVV_dB(self, **kwargs):
-        """Return VV backscattering coefficient in dB. Any parameter can be added to slice the results (e.g. frequency=37e9). See xarray slicing with sel method (to document)"""
+        """Return VV backscattering coefficient in dB. Any parameter can be added to slice the results (e.g. frequency=37e9).
+         See xarray slicing with sel method (to document)"""
         return dB(self.sigmaVV(**kwargs))
 
     def sigmaHH(self, **kwargs):
-        """Return HH backscattering coefficient. Any parameter can be added to slice the results (e.g. frequency=37e9). See xarray slicing with sel method (to document)"""
+        """Return HH backscattering coefficient. Any parameter can be added to slice the results (e.g. frequency=37e9).
+         See xarray slicing with sel method (to document)"""
         return self.sigma(polarization_inc='H', polarization='H', **kwargs)
 
-
     def sigmaHH_dB(self, **kwargs):
-        """Return HH backscattering coefficient in dB. Any parameter can be added to slice the results (e.g. frequency=37e9). See xarray slicing with sel method (to document)"""
+        """Return HH backscattering coefficient in dB. Any parameter can be added to slice the results (e.g. frequency=37e9).
+         See xarray slicing with sel method (to document)"""
         return dB(self.sigmaHH(**kwargs))
 
     def sigmaHV(self, **kwargs):
-        """Return HV backscattering coefficient. Any parameter can be added to slice the results (e.g. frequency=37e9). See xarray slicing with sel method (to document)"""
+        """Return HV backscattering coefficient. Any parameter can be added to slice the results (e.g. frequency=37e9).
+         See xarray slicing with sel method (to document)"""
         return self.sigma(polarization_inc='H', polarization='V', **kwargs)
 
     def sigmaHV_dB(self, **kwargs):
-        """Return HV backscattering coefficient in dB. Any parameter can be added to slice the results (e.g. frequency=37e9). See xarray slicing with sel method (to document)"""
+        """Return HV backscattering coefficient in dB. Any parameter can be added to slice the results (e.g. frequency=37e9).
+         See xarray slicing with sel method (to document)"""
         return dB(self.sigmaHV(**kwargs))
 
     def sigmaVH(self, **kwargs):
-        """Return VH backscattering coefficient. Any parameter can be added to slice the results (e.g. frequency=37e9). See xarray slicing with sel method (to document)"""
+        """Return VH backscattering coefficient. Any parameter can be added to slice the results (e.g. frequency=37e9).
+         See xarray slicing with sel method (to document)"""
         return self.sigma(polarization_inc='V', polarization='H', **kwargs)
 
     def sigmaVH_dB(self, **kwargs):
-        """Return VH backscattering coefficient in dB. Any parameter can be added to slice the results (e.g. frequency=37e9). See xarray slicing with sel method (to document)"""
+        """Return VH backscattering coefficient in dB. Any parameter can be added to slice the results (e.g. frequency=37e9).
+         See xarray slicing with sel method (to document)"""
         return dB(self.sigmaVH(**kwargs))
 
-    def save(self, filename):
-        """save a result to disk. Under the hood, this is a netCDF file produced by xarray (http://xarray.pydata.org/en/stable/io.html)."""
-        self.data.to_netcdf(filename)
-    #def groupby(self, variable):
+
+    # def groupby(self, variable):
     #    """iterated over a given variable. Variable is typically frequency, theta, polarization or snowpack"""
     #
     #    return ResultGroup(self.data.groupby(variable))
@@ -209,19 +263,26 @@ def save(self, filename):
 
 
 def concat_results(result_list, coord):
-    """Concatenate several results from :py:meth:`smrt.core.model.Model.run` (of type :py:class:`Result`) into a single result (of type :py:class:`Result`). This extends
-    the number of dimension in the xarray hold by the instance. The new dimension is specified with coord
+    """Concatenate several results from :py:meth:`smrt.core.model.Model.run` (of type :py:class:`Result`) into a single result
+    (of type :py:class:`Result`). This extends the number of dimension in the xarray hold by the instance. The new dimension
+    is specified with coord
 
     :param result_list: list of results returned by :py:meth:`smrt.core.model.Model.run` or other functions.
-    :param coord: a tuple (dimension_name, dimension_values) for the new dimension. Dimension_values must be a sequence or array with the same length as result_list.
+    :param coord: a tuple (dimension_name, dimension_values) for the new dimension. Dimension_values must be a sequence or
+    array with the same length as result_list.
 
     :returns: :py:class:`Result` instance
 
     """
 
     dim_name, dim_value = coord
 
-    return Result(xr.concat([result.data for result in result_list], pd.Index(dim_value, name=dim_name)))
+    ResultClass = type(result_list[0])
+
+    if not all([type(result) == ResultClass for result in result_list]):
+        raise SMRTError("The results are not all of the same type")
+
+    return ResultClass(xr.concat([result.data for result in result_list], pd.Index(dim_value, name=dim_name)))
 
 
 def _strongsqueeze(x):

smrt/core/test_result.py

@@ -7,12 +7,16 @@
 
 
 # Tests written in response to -ve intensity bug in result.py
-res_example = result.Result([[[[4.01445680e-03, 3.77746658e-03, 0.00000000e+00]],
+res_example = result.ActiveResult([[[[4.01445680e-03, 3.77746658e-03, 0.00000000e+00]],
                     [[3.83889082e-03, 3.85904771e-03, 0.00000000e+00]],
                     [[2.76453599e-20, -2.73266027e-20, 0.00000000e+00]]]],
                     coords = [('theta', [35]), ('polarization', ['V','H','U']),
                     ('theta_inc', [35]), ('polarization_inc', ['V','H','U'])])
 
+def test_methods():
+    assert hasattr(res_example, "sigma")
+    assert not hasattr(res_example, "Tb")
+
 def test_positive_sigmaVV():
     ok_(res_example.sigmaVV()>0)
 
@@ -37,3 +41,4 @@ def test_sigmaHV_dB():
 
 def test_sigmaVH_dB():
     np.testing.assert_allclose(res_example.sigmaVH_dB(), -14.0321985560285)
+

smrt/rtsolver/dort.py

@@ -26,7 +26,7 @@
 
 # local import
 from ..core.error import SMRTError
-from ..core.result import Result
+from ..core.result import make_result
 from ..core import lib
 from smrt.core.lib import smrt_matrix, isnull
 # Lazy import: from smrt.interface.coherent_flat import process_coherent_layers 
@@ -135,10 +135,11 @@ def solve(self, snowpack, emmodels, sensor, atmosphere=None):
             intensity = intensity.reshape(list(intensity.shape[:-1])+[intensity.shape[-1]//npol, npol])
         #  describe the results list of (dimension name, dimension array of value)
         coords = [('theta', sensor.theta_deg), ('polarization', pola)]
+
         if sensor.mode == 'A':
             coords = [('theta_inc', sensor.theta_inc_deg), ('polarization_inc', pola)] + coords
 
-        return Result(intensity, coords)
+        return make_result(sensor.mode, intensity, coords)
 
     def dort(self, m_max=0, special_return=False):
         # not to be called by the user
@@ -368,7 +369,7 @@ def dort_modem_banded(self, m, streams, eigenvalue_solver, interfaces, intensity
 
             # fill the vector
             if m == 0 and self.temperature is not None and self.temperature[l] > 0:
-                if Rtop_l is 0:
+                if isnull(Rtop_l):
                     b[il_topl:il_topl+nsl_npol, :] -= self.temperature[l]  # a mettre en (l)
                 else:
                     b[il_topl:il_topl+nsl_npol, :] -= ((1.0 - muleye(Rtop_l)) * self.temperature[l])[:, np.newaxis]  # a mettre en (l)