Jump to content
Fórum Script Brasil
  • 0

Erro ao executar CMOD5.N com python


nina_99

Question

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)
 
Edited by nina_99
Link to comment
Share on other sites

7 answers to this question

Recommended Posts

  • 0

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]

Link to comment
Share on other sites

  • 0

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]

 

Link to comment
Share on other sites

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.

Guest
Answer this question...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.



  • Forum Statistics

    • Total Topics
      152.1k
    • Total Posts
      652k
×
×
  • Create New...