Hello, i’m having trouble running the function cmod5n_inverse, i’ve inserted the parameters phi(wind direction-azimuth angle), and the incidence angle and sigma0 extracted from SNAP from a Sentinel-1 image, using linear regression.The function is supposed to return V, the wind velocity in meters per second, however the results seem to repeat themselves.Like this:
Initial guess V: [10. 10. 10. 10. 10.]
Iteration 1, V: [20. 20. 20. 20. 20.], step: 5.0, sigma0_calc: [0.00342216 0.00331692 0.00299704 0.00234962 0.00433662]
Iteration 2, V: [25. 25. 25. 25. 25.], step: 2.5, sigma0_calc: [0.02235289 0.02171421 0.01968306 0.01462111 0.02769178]
Iteration 3, V: [27.5 27.5 27.5 27.5 27.5], step: 1.25, sigma0_calc: [0.03657467 0.03551247 0.0321561 0.0236992 0.04577651]
Iteration 4, V: [28.75 28.75 28.75 28.75 28.75], step: 0.625, sigma0_calc: [0.0431854 0.04192145 0.0379445 0.0279204 0.05434188]
Final V after iterations: [28.75 28.75 28.75 28.75 28.75]
Wind Speed (m/s): [28.75 28.75 28.75 28.75 28.75]
Can someone inform me what is the possible error?
Here is the following code used to obtain the results:
!pip install cartopy
!pip install netCDF4
!pip install nansat
!pip install geopandas
import numpy as np
import warnings
warnings.simplefilter("ignore", RuntimeWarning)
def cmod5n_forward(v, phi, theta):
'''! ---------
! cmod5n_forward(v, phi, theta)
! inputs:
! v in [m/s] wind velocity (always >= 0)
! phi in [deg] angle between azimuth and wind direction
! (= D - AZM)
! theta in [deg] incidence angle
! output:
! CMOD5_N NORMALIZED BACKSCATTER (LINEAR)
!
! All inputs must be Numpy arrays of equal sizes
!
! A. STOFFELEN MAY 1991 ECMWF CMOD4
! A. STOFFELEN, S. DE HAAN DEC 2001 KNMI CMOD5 PROTOTYPE
! H. HERSBACH JUNE 2002 ECMWF COMPLETE REVISION
! J. de Kloe JULI 2003 KNMI, rewritten in fortan90
! A. Verhoef JAN 2008 KNMI, CMOD5 for neutral winds
! K.F.Dagestad OCT 2011 NERSC, Vectorized Python version
!---------------------------------------------------------------------
'''
from numpy import cos, exp, tanh, array
DTOR = 57.29577951
THETM = 40.
THETHR = 25.
ZPOW = 1.6
# NB: 0 added as first element below, to avoid switching from 1-indexing to 0-indexing
C = [0, -0.6878, -0.7957, 0.3380, -0.1728, 0.0000, 0.0040, 0.1103, 0.0159,
6.7329, 2.7713, -2.2885, 0.4971, -0.7250, 0.0450,
0.0066, 0.3222, 0.0120, 22.7000, 2.0813, 3.0000, 8.3659,
-3.3428, 1.3236, 6.2437, 2.3893, 0.3249, 4.1590, 1.6930]
Y0 = C[19]
PN = C[20]
A = C[19]-(C[19]-1)/C[20]
B = 1./(C[20]*(C[19]-1.)**(3-1))
# ! ANGLES
FI=phi/DTOR
CSFI = cos(FI)
CS2FI= 2.00 * CSFI * CSFI - 1.00
X = (theta - THETM) / THETHR
XX = X*X
# ! B0: FUNCTION OF WIND SPEED AND INCIDENCE ANGLE
A0 =C[ 1]+C[ 2]*X+C[ 3]*XX+C[ 4]*X*XX
A1 =C[ 5]+C[ 6]*X
A2 =C[ 7]+C[ 8]*X
GAM=C[ 9]+C[10]*X+C[11]*XX
S0 =C[12]+C[13]*X
# V is missing! Using V=v as substitute, this is apparently correct
V=v
S = A2*V
S_vec = S.copy()
SlS0 = [S_vec<S0]
S_vec[SlS0]=S0[SlS0]
A3=1./(1.+exp(-S_vec))
SlS0 = (S<S0)
A3[SlS0]=A3[SlS0]*(S[SlS0]/S0[SlS0])**( S0[SlS0]*(1.- A3[SlS0]))
#A3=A3*(S/S0)**( S0*(1.- A3))
B0=(A3**GAM)*10.**(A0+A1*V)
# ! B1: FUNCTION OF WIND SPEED AND INCIDENCE ANGLE
B1 = C[15]*V*(0.5+X-tanh(4.*(X+C[16]+C[17]*V)))
B1 = C[14]*(1.+X)- B1
B1 = B1/(exp( 0.34*(V-C[18]) )+1.)
# ! B2: FUNCTION OF WIND SPEED AND INCIDENCE ANGLE
V0 = C[21] + C[22]*X + C[23]*XX
D1 = C[24] + C[25]*X + C[26]*XX
D2 = C[27] + C[28]*X
V2 = (V/V0+1.)
V2ltY0 = V2<Y0
V2[V2ltY0] = A+B*(V2[V2ltY0]-1.)**PN
B2 = (-D1+D2*V2)*exp(-V2)
# ! CMOD5_N: COMBINE THE THREE FOURIER TERMS
CMOD5_N = B0*(1.0+B1*CSFI+B2*CS2FI)**ZPOW
return CMOD5_N
def cmod5n_inverse(sigma0_obs, phi, incidence, iterations=10):
from numpy import ones, array
# First guess wind speed
V = array([10.]) * ones(sigma0_obs.shape)
step = 10.
# Debugging prints
print(f"Initial guess V: {V}")
# Iterating until error is smaller than threshold
for iterno in range(1, iterations):
sigma0_calc = cmod5n_forward(V, phi, incidence)
ind = sigma0_calc - sigma0_obs > 0
V = V + step
V[ind] = V[ind] - 2 * step
step = step / 2
# Debugging prints
print(f"Iteration {iterno}, V: {V}, step: {step}, sigma0_calc: {sigma0_calc}")
# Debugging prints
print(f"Final V after iterations: {V}")
return V
# Data
sigma0_obs = [4.32674215, 3.72549955, 3.25660927, 2.93046451, 4.58418151]
incidence = [64.1597756, 64.8628612, 67.0776741, 72.8808313, 57.7738178]
phi = [-261.571387, -262.019219, -264.074833, -267.473085, -270.674255]
# Convert lists to numpy arrays
sigma0_obs = np.array(sigma0_obs)
incidence = np.array(incidence)
phi = np.array(phi)
# Function call to cmod5n_inverse to get wind speed (V)
V = cmod5n_inverse(sigma0_obs, phi, incidence, iterations=5)
# Print the resulting wind speed (V)
print("Wind Speed (m/s):", V)