Example 2: Forward modeling of RT model combination (surface + canopy part)

1. Narrative (Examples of different model combination (surface + canopy part))

This notebook demonstrates model combinations of

Dubois95 + SSRT
Oh04 + WaterCloud

A short overview of the required input parameters can be found here.

2. Requirements

Installation of SenSE

3. Required imports

import numpy as np
from sense.util import f2lam
from sense.model import RTModel
from sense.soil import Soil
from sense.canopy import OneLayer
import matplotlib.pyplot as plt

4. Dubois95+SSRT for different incidence angles

4.1 Input parameters for the RT model combination

# array with different incidence angles in radians
theta_deg = np.arange(0.,70.)
theta = np.deg2rad(theta_deg) # incidence angle [radians]

# soil model parameters
f  = 13.  # frequency [GHz]
lam = f2lam(f)  # wavelength [m]
s = 0.15/100.  # surface roughness [m]
eps = 15. - 4.0j # dielectric constant

# canopy parameters for short alfalfa
d = 0.17 # vegetation height [m]
tau = 2.5 # optical depths
ke = tau/d # extinction coefficient [m⁻¹]
omega = 0.27 # scattering albedo
ks=omega*ke

4.2 Forward modeling of RT model combination Dubois95 and SSRT (rayleigh scattering)

First model run with model parameters for small alfalfa

# Define models and polarization to be used
models = {'surface' : 'Dubois95', 'canopy' : 'turbid_rayleigh'} # For SSRT two scattering models are implemented 'turbid_rayleigh' and 'turbid_isotropic'
# models = {'surface' : 'Oh92', 'canopy' : 'turbid_isotropic'} # alternative with different surface and canopy models
pol='vv' # target polarization

# Soil model initialization based on previously defined input parameters
S = Soil(f=f, s=s, eps=eps)

# Canopy model initialization based on previously defined input parameters
C = OneLayer(ke_h=ke, ke_v=ke, d=d, ks_v=ks, ks_h=ks, canopy=models['canopy'])

# Combined Model initialization
RT = RTModel(theta=theta, models=models, surface=S, canopy=C, freq=f)

# Run RT model
RT.sigma0()
back_short = RT.stot[pol]

/home/docs/checkouts/readthedocs.org/user_builds/sense-community-sar-scattering-model/envs/dev/lib/python3.10/site-packages/sense/surface/dubois95.py:42: RuntimeWarning: divide by zero encountered in divide
  b = 10**-2.35 * (np.cos(self.theta) ** 3) / (np.sin(self.theta) ** 3)
/home/docs/checkouts/readthedocs.org/user_builds/sense-community-sar-scattering-model/envs/dev/lib/python3.10/site-packages/sense/surface/dubois95.py:45: RuntimeWarning: invalid value encountered in multiply
  return b * c * d
/home/docs/checkouts/readthedocs.org/user_builds/sense-community-sar-scattering-model/envs/dev/lib/python3.10/site-packages/sense/surface/dubois95.py:35: RuntimeWarning: divide by zero encountered in divide
  a = 10**-2.75 * (np.cos(self.theta) ** 1.5) / (np.sin(self.theta) ** 5)
/home/docs/checkouts/readthedocs.org/user_builds/sense-community-sar-scattering-model/envs/dev/lib/python3.10/site-packages/sense/surface/dubois95.py:38: RuntimeWarning: invalid value encountered in multiply
  return a * c * d

Second model run with model parameters (canopy) for tall alfalfa

# redefine canopy parameters for tall alfalfa
d = 0.55 # vegetation height [m]
tau = 0.45 # optical depths
ke = tau/d # extinction coefficient [m⁻¹]
omega = 0.175 # scattering albedo
ks=omega*ke

# initialize and run of the model combination
S = Soil(f=f, s=s, eps=eps)
C = OneLayer(ke_h=ke, ke_v=ke, d=d, ks_v=ks, ks_h=ks, canopy=models['canopy'])
RT = RTModel(theta=theta, models=models, surface=S, canopy=C, freq=f)
RT.sigma0()

4.3 Visualization

# plot backscatter changes for short and high alfalfa based on the incidence angle
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(theta_deg, 10.*np.log10(back_short), label='short alfalfa', color='b')
ax.plot(theta_deg, 10.*np.log10(RT.stot[pol]), label='tall alfalfa', color='r')

ax.legend()

ax.grid()
ax.set_xlabel('Incidence angle [°]')
ax.set_ylabel('Sigma vv [dB]')
ax.set_xlim(0.,70.)
ax.set_ylim(-16.,6.)

plt.show()

../_images/47b7a0f4b8379633f1da801673d87d43c3e4709b6df260c09d4c29a5b6570546.png

5. Oh04+WaterCloud for different incidence angles

5.1 Input parameters for the RT model combination

# Parameters for both models
theta = np.deg2rad(40) # incidence angle [radians]
f  = 5.3 # Frequency [GHz]

# soil parameter
s = 0.3/100. # surface roughness [m]
mv = np.linspace(0.01,0.4) # soil moisture [m³/m³]
# calculations
lam = f2lam(f)  # wavelength [m]
k = 2.*np.pi/lam # radar wave number
ks = k * s

# canopy parameters
A_vv = 0.0950 # empirical parameter - need to be optimized for each test site
B_vv = 0.5513 # empirical parameter - need to be optimized for each test site
A_hh = 0 # empirical parameter - need to be optimized for each test site
B_hh = 0 # empirical parameter - need to be optimized for each test site
A_hv = 0 # empirical parameter - need to be optimized for each test site
B_hv = 0 # empirical parameter - need to be optimized for each test site

V1 = 0.02 # parameter describing the canopy (todo: find realistic values)
V2 = np.linspace(0.5,5) # # parameter describing the canopy (todo: find realistic values)

5.2 Forward modeling of RT model combination Oh04 and WaterCloud

# Define models and polarization to be used
models = {'surface' : 'Oh04', 'canopy' : 'water_cloud'}
pol='vv' # target polarization

# Soil model initialization based on previously defined input parameters
S = Soil(f=f, s=s, mv=mv)

# Canopy model initialization based on previously defined input parameters
C = OneLayer(A_vv=A_vv, B_vv=B_vv, A_hh=A_hh, B_hh=B_hh, A_hv=A_hv, B_hv=B_hv, V1=V1, V2=V2, d=1, ks_v=1, ke_v=1, ke_h=1, canopy=models['canopy'])

# Combined Model initialization
RT = RTModel(theta=theta, models=models, surface=S, canopy=C, freq=f)

# Run RT model
RT.sigma0()
backscatter_vv = RT.stot[pol] # total modeled backscatter
backscatter_vv_g = RT.s0g[pol] # modeled ground contribution
backscatter_vv_c = RT.s0c[pol] # modeled canopy contribution

WARNING: Permittivity cannot be calculated due to missing soil texture!

5.3 Visualization

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(mv, 10.*np.log10(backscatter_vv), label='total backscatter', color='b')
ax.plot(mv, 10.*np.log10(backscatter_vv_g), label='ground contribution', color='r')
ax.plot(mv, 10.*np.log10(backscatter_vv_c), label='canopy contribution', color='g')

ax.legend()

ax.grid()
ax.set_xlabel('Soil moisture [m³/m³]')
ax.set_ylabel('Sigma nought [dB]')

ax.set_ylim(-40.,-10.)

plt.show()

../_images/4bcbf0e8ea1eba560f8df09626dff2ab2cfd888e662539a55dc11c7722fd2202.png