smrt.rtsolver package

Contains different solvers of the radiative transfer equation. Based on the electromagnetic properties of each layer computed by the EM model, these RT solvers compute the emission and propagation of energy in the medium up to the surface (the atmosphere is usually dealt with independently in dedicated modules in smrt.atmosphere).

The solvers differ by the approximations and numerical methods. dort is currently the most accurate and recommended in most cases unless the computation time is a constraint.

Selection of the solver is done with the smrt.core.model.make_model function.

Here are some recommendations to chose the appropriate solver:

SMRT Solvers Overview

Solver Name	Description	Use Case	Performance
dort	The default original solver in SMRT for passive microwave (brightness temperature) and radar (sigma0). General and robust.	General applications for both passive microwave and radar.	Robust but slow
iterative_first_order	Radar-only solver using an iterative method up to first order. Fast and provides contributions of interaction mechanisms.	Radar studies, especially to analyze interaction mechanisms.	Fast
multifresnel_thermalemission	Passive microwave-only solver assuming no scattering in the layers and flat interfaces.	Radiometry at low frequency when scattering is negligible.	Very fast
nadir_lrm_altimetry	Specialized solver for altimeter in Low Rate Mode. Solve the first order interaction.	Altimetry in Low Rate Mode (LRM).

Note: For Developers

To experiment with DORT, it is recommended to copy the file dort.py to e.g. dort_mytest.py so it is immediately available through smrt.core.model.make_model.

To develop a new solver that will be accessible by the smrt.core.model.make_model function, add a file in this directory. Refer to dort.py as an example. Only the solve method needs to be implemented. It must return a smrt.core.result.Result instance with the results. Contact the core developers for more details.

smrt.rtsolver.dort module

Provides the Discrete Ordinate and Eigenvalue Solver as a multi-stream solver of the radiative transfer model.

This solver is precise but less efficient than 2 or 6 flux solvers. Different flavours of DORT (or DISORT) exist in the literature depending on the mode (passive or active), on the density of the medium (sparse media have trivial inter-layer boundary conditions), on the way the streams are connected between the layers and on the way the phase function is prescribed. The actual version is a blend between Picard et al. 2004 (active mode for sparse media) and DMRT-ML (Picard et al. 2013) which works in passive mode only for snow. The DISORT often used in optics (Stamnes et al. 1988) only works for sparse medium and uses a development of the phase function in Legendre polynomials on theta. The version used in DMRT-QMS (L. Tsang’s group) is similar to the present implementation except it uses spline interpolation to connect constant-angle streams between the layers although we use direct connection by varying the angle according to Snell’s law. A practical consequence is that the number of streams vary (due to internal reflection) and the value n_max_stream only applies in the most refringent layer. The number of outgoing streams in the air is usually smaller, sometimes twice smaller (depends on the density profile). It is important not to set too low a value for n_max_stream. E.g. 32 is usually fine, 64 or 128 are better but simulations will be much slower.

Note

The DORT solver is very robust in passive mode but may raise exception in active mode due to a matrix diagonalisation problem. The exception provides detailed information on how to address this issue. Two new diagonalisation approaches were added in January 2024. They are activated by setting the diagonalization_method optional argument (see smrt.core.make_model). The first method (diagonalization_method=’shur’) replaces the scipy.linalg.eig function by a shur decomposition followed by a diagonalisation of the shur matrix. While scipy.linalg.eig performs such a shur decomposition internally in any case, it seems that explicitly calling the shur decomposition beforehand improves the stability. Nevertheless to really solve the problem, the second method (diagonalization_method=’shur_forcedtriu’) consists in removing the 2x2 and 3x3 blocks from the shur matrix, i.e. forcing the shur matrix to be upper triangular (triu in numpy jargon = zeroing the lower part of this matrix). This problem is due to the structure of the matrix to be diagonalized and the formulation of the DORT method in the polarimetric configuration. The eigenvalues come by triplets and can be very close to each other for the three H, V, U Stokes components when scattering is becoming small (or equiv. the azimuth mode ‘m’ is large). As a consequence of the Gershgorin theorem, this results in slightly complex eigenvalues (i.e. eigenvalues with very small imaginary part) that comes from 2x2 or 3x3 blocks in the shur decomposition. This would not be a problem if the eigenvectors were correctly estimated, but this is not the case. It is indeed difficult to find the correct orientation of eigenvectors associated to very close eigenvalues. To overcome the problem, the solution is to remove the 2x2 and 3x3 blocks. In principle, it would be safer to check that these blocks are nearly diagonal but this is not done in the current implementation. The user is responsible to switch between the options until it works. After sufficient successful reports by users are received the last method (forcedtriu) will certainly be the default.

class DORT(n_max_stream=32, m_max=2, stream_mode='most_refringent', phase_normalization='auto', error_handling='exception', process_coherent_layers=False, prune_deep_snowpack=None, diagonalization_method='eig')

Bases: object

Implements the Discrete Ordinate and Eigenvalue Solver.

Parameters:

• n_max_stream – number of stream in the most refringent layer.

• m_max – number of mode (azimuth).

• phase_normalization – the integral of the phase matrix should in principe be equal to the scattering coefficient. However, some emmodels do not respect this strictly. In general a small difference is due to numerical rounding and is acceptable, but a large difference rather indicates either a bug in the emmodel or input parameters that breaks the assumption of the emmodel. The most typical case is when the grain size is too big compared to wavelength for emmodels that rely on Rayleigh assumption. If this argument is to True, the phase matrix is normalized to be coherent with the scattering coefficient, but only when the difference is moderate (0.7 to 1.3). If set to “force” the normalization is always performed. This option is dangerous because it may hide bugs or unappropriate input parameters (typically too big grains). If set to False, no normalization is performed. If set to “auto” the normalization is performed except for emmodels not respecting the reciprocity princple (which the normalization relies on).

• stream_mode – TODO: add documentation.

• error_handling – If set to “exception” (the default), raise an exception in case of error, stopping the code. If set to “nan”, return a nan, so the calculation can continue, but the result is of course unusuable and the error message is not accessible. This is only recommended for long simulations that sometimes produce an error.

• process_coherent_layers – Adapt the layers thiner than the wavelegnth using the MEMLS method. The radiative transfer theory is inadequate layers thiner than the wavelength and using DORT with thin layers is generally not recommended. In some parcticular cases (such as ice lenses) where the thin layer is isolated between large layers, it is possible to replace the thin layer by an equivalent reflective interface. This neglects scattering in the thin layer, which is acceptable in most case, because the layer is thin. To use this option and more generally to investigate ice lenses, it is recommended to read MEMLS documentation on this topic.

• prune_deep_snowpack – this value is the optical depth from which the layers are discarded in the calculation. It is to be use to accelerate the calculations for deep snowpacks or at high frequencies when the contribution of the lowest layers is neglegible. The optical depth is a good criteria to determine this limit. A value of about 6 is recommended. Use with care, especially values lower than 6.

solve(snowpack, emmodels, sensor, atmosphere=None)

Solve the radiative transfer equation for a given snowpack, emmodels and sensor configuration.

Parameters:

• snowpack – Snowpack object.

• emmodels – List of electromagnetic models.

• sensor – Sensor object.

• atmosphere – Optional atmosphere object.

Returns:

Computed result.

Return type:

Result object

smrt.rtsolver.dort_nonormalization module

Provides a Discrete Ordinate and Eigenvalue Solver as a multi-stream solver of the radiative transfer model.

Solvers are precise but less efficient than 2 or 6 flux solvers. Different flavours of DORT (or DISORT) exist depending on the mode (passive or active), on the density of the medium (sparse media have trivial inter-layer boundary conditions), on the way the streams are connected between the layers and on the way the phase function is prescribed. The actual version is a blend between Picard et al. 2004 (active mode for sparse media) and DMRT-ML (Picard et al. 2013) which works in passive mode only for snow. The DISORT often used in optics (Stamnes et al. 1988) works only for sparse medium and uses a development of the phase function in Legendre polynomia on theta. The version used in DMRT-QMS (L. Tsang’s group) is similar to the present implementation except it uses spline interpolation to connect constant-angle streams between the layers although we use direct connection by varying the angle according to Snell’s law. A practical consequence is that the number of streams vary (due to internal reflection) and the value n_max_stream only applies in the most refringent layer. The number of outgoing streams in the air is usually smaller, sometimes twice smaller (depends on the density profile). It is important not to set too low a value for n_max_stream. E.g. 32 is usually fine, 64 or 128 are better but simulations will be much slower.

class DORT(n_max_stream=32, m_max=2, stream_mode='most_refringent')

Bases: object

Implements the Discrete Ordinate and Eigenvalue Solver.

Parameters:

• n_max_stream (int) – Number of streams in the most refringent layer.

• m_max (int) – Number of mode (azimuth).

• stream_mode (str) – Stream mode, default is “most_refringent”.

solve(snowpack, emmodels, sensor, atmosphere=None)

Solves the radiative transfer equation for a given snowpack, emmodels and sensor configuration.

Parameters:

• snowpack – Snowpack object.

• emmodels – List of electromagnetic models.

• sensor – Sensor object.

• atmosphere – Optional atmosphere object.

Returns:

Computed result.

Return type:

Result

smrt.rtsolver.nadir_lrm_altimetry module

Computes waveforms as measured by Low Rate Mode altimeters (e.g. ENVISAT) for the given snowpack and sensor (or complex terrain soon). The implementation is based on Adams and Brown 1998 and Lacroix et al. 2008. Both models differ in the specific choices for the scattering and backscatter of the interface, but are similar in the way the waveform is calculated, which concerns the solver here.

Approximations:

• Backscatter is computed assuming only first order scattering. The propagation is then simply governed by extinction.

• Near nadir / small angle approximation: to compute delay, the paths in the snow are along the z-axis. Off-nadir delay is neglected. This error is likely to be small (except for very deep penetration).

Note

With this RT solver, if using Geometrical Optics for rough surface/interface modeling, it is strongly advised to use smrt.interface.geometrical_optics_backscatter() instead of smrt.interface.geometrical_optics() for the reason explained in the documentation of those modules.

class NadirLRMAltimetry(waveform_model=None, oversampling=10, return_oversampled=False, skip_pfs_convolution=False, return_contributions=False, theta_inc_sampling=8, return_theta_inc_sampling=False, error_handling='exception')

Bases: object

Implements the Nadir LRM Mode Altimetry RT solver.

Parameters:

• oversampling – integer number defining the number of subgates used for the computation in each altimeter gate. This is equivalent to multiply the bandwidth by this number. It is used to perform more accurate computation.

• return_oversampled – by default the backscatter is returned for each gate. If set to True, the oversampled waveform is returned instead. See the ‘oversampling’ argument.

• return_contributions – return volume, surface and interface backscatter contributions in addition to the total backscatter.

• skip_pfs_convolution – return the vertical backscatter without the convolution by the PFS, if set to True.

• theta_inc_sampling – number of subdivisions used to calculate the incidence angular variations of surface and inteface backscatter (the higher the better but the more computationnaly expensive). Note that the subdivisions are irregular in incidence angle but correspond to annulii of equi-duration. This number must be a true divider of the number of gates.

• return_theta_inc_sampling – return the backscatter at the different angles

• error_handling – If set to “exception” (the default), raise an exception in case of error, stopping the code. If set to “nan”, return a nan, so the calculation can continue, but the result is of course unusuable and the error message is not accessible. This is only recommended for long simulations that sometimes produce an error.

solve(snowpack, emmodels, sensor, atmosphere=None)

Solves the radiative transfer equation for a given snowpack, emmodels and sensor configuration.

Parameters:

• snowpack – Snowpack object.

• emmodels – List of electromagnetic models.

• sensor – Sensor object.

• atmosphere – Optional atmosphere object.

Returns:

Computed result.

Return type:

AltimetryResult

gate_depth(eps=None)

return gate depth that cover the snowpack for a regular time sampling

vertical_scattering_distribution(return_contributions, mu_i=1.0)

Compute the vertical backscattering distribution due to “grain” or volume scattering (symbol pvg in Eq 9 in Lacroix 2008) and “interfaces” or ‘surface’ scattering (symbol pvl in Eq 9 in Lacroix 2008)

Parameters:

mu – cosine of the incidence angles. Only the dependence on the surface scattering depend on mu_i

Returns:

Backscattering distribution.

Return type:

ndarray

smrt.rtsolver.waveform_model module

class Brown1977(sensor, numerical_convolution=False)

Bases: WaveformModel

Implements the Antenna Gain formulation used by Brown 1977.

The formula is exp(2/gamma * sin(theta)**2) for the perfect nadir case, but is also available with off-nadir angles.

Parameters:

• sensor – Sensor object.

• numerical_convolution (bool) – Whether to use numerical convolution.

PFS_PTR_PDF(tau, sigma_surface=0, surface_slope=0)

Computes the convolution of the PFS and PTR.

Parameters:

• tau – Time to which to compute the PFSxPTR.

• sigma_surface (float) – RMS height of the surface topography (meter).

• surface_slope (float) – Surface slope.

Returns:

Convoluted PFS and PTR.

Return type:

ndarray

class Newkrik1992(sensor)

Bases: WaveformModel

Implements the Antenna Gain formulation proposed by Newkrik and Brown, 1992.

Compared to the classical Brown 1977, it takes into account the asymmetry of the antenna pattern in the co and cross-track direction.

smrt.rtsolver.iterative_first module