Error ruinning CMOD5.N function with python

29 views Asked by At

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





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)
0

There are 0 answers