Erro ao executar CMOD5.N com python

[nina_99]

Olá, estou tendo problemas para executar a função cmod5n_inverse, inseri os parâmetros phi(direção do vento-ângulo azimutal), e o ângulo de incidência e sigma0 extraídos do SNAP de uma imagem do Sentinel-1, usando regressão linear. A função é deveria retornar V, a velocidade do vento em metros por segundo, porém os resultados parecem se repetir. Assim:

Suposição inicial V: [10. 10. 10. 10. 10.]
Iteração 1, V: [20. 20. 20. 20. 20.], passo: 5.0, sigma0_calc: [0.00342216 0.00331692 0.00299704 0.00234962 0.00433662]
Iteração 2, V: [25. 25. 25. 25. 25.], passo: 2.5, sigma0_calc: [0.02235289 0.02171421 0.01968306 0.01462111 0.02769178]
Iteração 3, V: [27,5 27,5 27,5 27,5 27,5], passo: 1,25, sigma0_calc: [0,03657467 0,03551247 0,0321561  0,0236992  0,04577651]
Iteração 4, V: [28,75 28,75 28,75 28,75 28,75], passo: 0,625, sigma0_calc: [0,0431854  0,04192145 0,0379445  0,0279204  0,05434188]
Final V após iterações: [28,75 28,75 28,75 28,75 28,75]
Velocidade do vento (m/s): [28,75 28,75 28,75 28,75 28,75]

Alguém pode me informar qual é o possível erro?
Aqui está o seguinte código usado para obter os resultados:

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]

 

sigma0_obs = np.array(sigma0_obs)

incidence=np.array(incidence)

phi=np.array(phi)

 

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

 

# Call cmod5n_inverse to get the wind speed (V) for the given observation

V = cmod5n_inverse(sigma0_obs,phi,incidence, iterations=5)

 

# Print the resulting wind speed (V)

print("Wind Speed (m/s):", V)

 

Editado Agosto 2, 2023 por nina_99

[ArteEN]

poderia colocar o codigo dentro das tags de codigo? assim aparecera organizado para copiar na formatação correta, mais confiavel

[nina_99]

Obrigada, adicionei as tags.

[ArteEN]

me referia a estas tags que formatam o codigo no forum

var = 1
var2 = [2, 3]

for i in var2:
  print(i + var)

 

[nina_99]

Pode me passar um link de como editar  o codigo neste formato?

[ArteEN]

ok, consegui executar e obter o mesmo resultado que você passou. agora qual seria o resultado esperado?

[nina_99]

Boa tarde, a velocidade deveria ser um so resultado, no exemplo acima esta repetindo o mesmo valor 5 vezes dentro do colchete.

Deveria sair desta forma:

Initial guess V: [10.]

Iteration 1, V: [20.], step: 5.0, sigma0_calc: [0.00436189]

Iteration 2, V: [25.], step: 2.5, sigma0_calc: [0.02770691]

Iteration 3, V: [27.5], step: 1.25, sigma0_calc: [0.0456733]

Iteration 4, V: [28.75], step: 0.625, sigma0_calc: [0.05417324]

Iteration 5, V: [29.375], step: 0.3125, sigma0_calc: [0.05813452]

Final V after iterations: [29.9609375] Wind Speed (m/s): [29.9609375]

[ArteEN]

veja a anotação dentro da função

Em 01/08/2023 em 21:49, nina_99 disse:

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

    !---------------------------------------------------------------------

       '''

diz que todas as entradas deve ser arrays de tamanho iguais, agora veja as entradas

Em 01/08/2023 em 21:49, nina_99 disse:

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]

são arrays de mesmo tamanho, o que esta acontecendo é que ele esta mostrando 5 resultados para 5 valores de dados, se você quizer apenas uma saida, deve então passar arrays com apenas um valor.

 

alterando os arrays para

sigma0_obs=[4.32674215]
incidence=[64.1597756]
phi=[-261.571387]

você terá a seguinte saida

Initial guess V: [10.]
Iteration 1, V: [20.], step: 5.0, sigma0_calc: [0.00342216]
Iteration 2, V: [25.], step: 2.5, sigma0_calc: [0.02235289]
Iteration 3, V: [27.5], step: 1.25, sigma0_calc: [0.03657467]
Iteration 4, V: [28.75], step: 0.625, sigma0_calc: [0.0431854]
Final V after iterations: [28.75]
Wind Speed (m/s): [28.75]