Error in running CMOD5.N with Python

ufcgeo · August 1, 2023, 11:14pm

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)

raavikant · September 7, 2023, 6:07am

Hi,

As per my knowledge, it seems like your code is trying to iteratively estimate the wind speed (V) using the cmod5n_inverse function based on observed sigma0 values, incidence angles, and azimuth angles (phi). However, you mentioned that the results are not as expected, and the values seem to be repeating.

You may try a few potential issues and suggestions to help you debug:

Initial Guess: Your initial guess for wind speed is set to 10 m/s for all observations. Depending on the problem, this may or may not be a suitable initial guess. You might want to experiment with different initial values.
Convergence Criteria: The loop in your cmod5n_inverse function iterates for a fixed number of times (10 iterations). In practice, you may want to implement a convergence criterion to stop the iteration when the change in V becomes very small, indicating convergence.
Step Size: The step size for adjusting V decreases by half in each iteration. If the change is too small or too large, it might lead to convergence issues. You can experiment with different step sizes or adapt the step size dynamically based on the convergence behavior.
Data Scaling: Ensure that the input data (sigma0_obs, incidence, phi) is correctly scaled and compatible with the assumptions made in the cmod5n_forward function. Any mismatch in units or assumptions could lead to convergence problems.
Debugging Output: Your code includes debugging print statements, which is a good practice. Review the output of these print statements to understand how the variables (V, step, sigma0_calc) change in each iteration and whether they behave as expected.
Check cmod5n_forward Function: Ensure that the cmod5n_forward function is implemented correctly and returns meaningful results. Check the formulas and constants used in this function to verify if they match the expected behavior.
Input Data: Double-check that the input data (sigma0_obs, incidence, phi) corresponds to the problem you are trying to solve. Any discrepancies in the input data can lead to unexpected results.
Initial Conditions: Make sure that the initial conditions and parameters for the problem are set correctly. For example, check if the constants (C) used in the cmod5n_forward function are appropriate for your specific application. You may also get assistance from an individual, having Python or mlops certification.

By addressing these points and experimenting with different initial guesses and convergence criteria, you should be able to debug and improve the performance of your wind speed estimation code.

Thanks

ufcgeo · December 14, 2023, 1:08pm

Thank you for the well thought out answer!