Skip to content

Commit

Permalink
Merge pull request #63 from vincelhx/branch_dev_vinc
Browse files Browse the repository at this point in the history 
added gmf noaa for RCM
  • Loading branch information
@agrouaze
agrouaze authored Oct 17, 2023
2 parents 8615d87 + 2327091 commit 7b685bb
Show file tree
Hide file tree
Showing 2 changed files with 172 additions and 94 deletions.
170 changes: 114 additions & 56 deletions src/xsarsea/windspeed/gmfs_impl.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,6 @@ def gmf_cmod5_generic(neutral=False):
6.2437, 2.3893, 0.3249, 4.159, 1.693])
name = 'gmf_cmod5n'


@GmfModel.register(name, wspd_range=[0.2, 50.], pol='VV', units='linear')
def gmf_cmod5(inc, wspd, phi):
zpow = 1.6
Expand Down Expand Up @@ -56,7 +55,8 @@ def gmf_cmod5(inc, wspd, phi):
b0 = (a3 ** gam) * 10. ** (a0 + a1 * wspd)

# B1 term
b1 = c[15] * wspd * (0.5 + x - np.tanh(4. * (x + c[16] + c[17] * wspd)))
b1 = c[15] * wspd * \
(0.5 + x - np.tanh(4. * (x + c[16] + c[17] * wspd)))
b1 = (c[14] * (1. + x) - b1) / (np.exp(0.34 * (wspd - c[18])) + 1.)

# B2 term
Expand All @@ -75,11 +75,13 @@ def gmf_cmod5(inc, wspd, phi):

return gmf_cmod5


# register gmfs gmf_cmod5 and gmf_cmod5n
gmf_cmod5_generic(neutral=False)
gmf_cmod5_generic(neutral=True)

@GmfModel.register( wspd_range=[0.2, 50.], pol='VV', units='linear')

@GmfModel.register(wspd_range=[0.2, 50.], pol='VV', units='linear')
def gmf_cmodifr2(inc_angle, wind_speed, wind_dir):

C = np.zeros(26)
Expand Down Expand Up @@ -141,8 +143,8 @@ def gmf_cmodifr2(inc_angle, wind_speed, wind_dir):
pt2 = 2 * tetanor * pt1 - pt0
# pt3 = 2 * tetanor * pt2 - pt1
b1 = C[8] + C[9] * pv1 \
+ (C[10] + C[11] * pv1) * pt1 \
+ (C[12] + C[13] * pv1) * pt2
+ (C[10] + C[11] * pv1) * pt1 \
+ (C[12] + C[13] * pv1) * pt2
tetamin = 18.0
tetamax = 58.0
tetanor = (2.0 * T - (tetamin + tetamax)) / (tetamax - tetamin)
Expand All @@ -158,12 +160,12 @@ def gmf_cmodifr2(inc_angle, wind_speed, wind_dir):
pt2 = 2 * tetanor * pt1 - pt0
# pt3 = 2 * tetanor * pt2 - pt1
result = (
C[14]
+ C[15] * pt1
+ C[16] * pt2
+ (C[17] + C[18] * pt1 + C[19] * pt2) * pv1
+ (C[20] + C[21] * pt1 + C[22] * pt2) * pv2
+ (C[23] + C[24] * pt1 + C[25] * pt2) * pv3
C[14]
+ C[15] * pt1
+ C[16] * pt2
+ (C[17] + C[18] * pt1 + C[19] * pt2) * pv1
+ (C[20] + C[21] * pt1 + C[22] * pt2) * pv2
+ (C[23] + C[24] * pt1 + C[25] * pt2) * pv3
)
b2 = result

Expand All @@ -176,8 +178,8 @@ def gmf_cmodifr2(inc_angle, wind_speed, wind_dir):
#
# gmf_dummy example
#
#@GmfModel.register(wspd_range=[3., 80.], pol='VH', units='linear')
#def gmf_dummy(incidence, speed, phi=None):
# @GmfModel.register(wspd_range=[3., 80.], pol='VH', units='linear')
# def gmf_dummy(incidence, speed, phi=None):
# a0 = 0.00013106836021008122
# a1 = -4.530598283705591e-06
# a2 = 4.429277425062766e-08
Expand All @@ -190,14 +192,13 @@ def gmf_cmodifr2(inc_angle, wind_speed, wind_dir):
# sig = a * speed ** b
#
# return sig



@GmfModel.register(wspd_range=[3., 80.],pol='VH', units='linear')
@GmfModel.register(wspd_range=[3., 80.], pol='VH', units='linear')
def gmf_rs2_v2(incidence, speed, phi=None):
"""
Radarsat-2 VH GMF : relation between sigma0, incidence and windspeed.
Parameters
----------
incidence: xarray.DataArray
Expand All @@ -209,49 +210,50 @@ def gmf_rs2_v2(incidence, speed, phi=None):
-------
sigma0: xarray.DataArray
linear sigma0
"""
#constants params

Z1_p = np.array([ 6.55519203e-06, 2.49753154e+00, -1.35734881e-02])
Z2_p = np.array([ 1.47342197e-04, -4.07334797e-06, 3.43593382e-08, 1.10188639e+00,
1.40782758e-02, -1.53748743e-04])
Final_p = np.array([-0.18675905, 24.48859492, 0.19185442, 25.38275738])
# constants params

#Z1
Z1_p = np.array([6.55519203e-06, 2.49753154e+00, -1.35734881e-02])
Z2_p = np.array([1.47342197e-04, -4.07334797e-06, 3.43593382e-08, 1.10188639e+00,
1.40782758e-02, -1.53748743e-04])
Final_p = np.array([-0.18675905, 24.48859492, 0.19185442, 25.38275738])

# Z1

a0_Z1 = Z1_p[0]
b0_Z1 = Z1_p[1]
b1_Z1 = Z1_p[2]
a_Z1 = a0_Z1
b_Z1 = b0_Z1 + b1_Z1*incidence
sig_Z1 = a_Z1*speed**(b_Z1)
#Z2

a_Z1 = a0_Z1
b_Z1 = b0_Z1 + b1_Z1*incidence
sig_Z1 = a_Z1*speed**(b_Z1)

# Z2
a0_Z2 = Z2_p[0]
a1_Z2 = Z2_p[1]
a2_Z2 = Z2_p[2]
b0_Z2 = Z2_p[3]
b1_Z2 = Z2_p[4]
b2_Z2 = Z2_p[5]
a_Z2 = a0_Z2 + a1_Z2*incidence + a2_Z2*incidence**2

a_Z2 = a0_Z2 + a1_Z2*incidence + a2_Z2*incidence**2
b_Z2 = b0_Z2 + b1_Z2*incidence + b2_Z2*incidence**2
sig_Z2 = a_Z2*speed**(b_Z2)

c0,c1 = Final_p[0],Final_p[1]
c2,c3 = Final_p[2],Final_p[3]
c0, c1 = Final_p[0], Final_p[1]
c2, c3 = Final_p[2], Final_p[3]
sigmoid_Z1 = 1 / (1 + np.exp(-c0*(speed-c1)))
sigmoid_Z2 = 1 / (1 + np.exp(-c2*(speed-c3)))
sig_Final = sig_Z1 * sigmoid_Z1 + sig_Z2 * sigmoid_Z2
sig_Final = sig_Z1 * sigmoid_Z1 + sig_Z2 * sigmoid_Z2
return sig_Final

@GmfModel.register(wspd_range=[3., 80.],pol='VH', units='linear')

@GmfModel.register(wspd_range=[3., 80.], pol='VH', units='linear')
def gmf_s1_v2(incidence, speed, phi=None):
"""
Sentinel-1 VH GMF : relation between sigma0, incidence and windspeed.
Parameters
----------
incidence: xarray.DataArray
Expand All @@ -264,40 +266,96 @@ def gmf_s1_v2(incidence, speed, phi=None):
sigma0: xarray.DataArray
linear sigma0
"""

#constants params

Z1_p = np.array([ 2.13755392e-06, 2.47395267e+00, -2.85775085e-03])

Z2_p = np.array([ 6.54058552e-05, -2.43845137e-06, 2.87698338e-08, 1.14509104e+00,
3.41828829e-02, -4.79715441e-04])
# constants params

Z1_p = np.array([2.13755392e-06, 2.47395267e+00, -2.85775085e-03])

Z2_p = np.array([6.54058552e-05, -2.43845137e-06, 2.87698338e-08, 1.14509104e+00,
3.41828829e-02, -4.79715441e-04])
Final_p = np.array([-0.23257086, 12.39717002, 0.21667263, 12.22862991])


#Z1
# Z1
a0_Z1 = Z1_p[0]
b0_Z1 = Z1_p[1]
b1_Z1 = Z1_p[2]

a_Z1 = a0_Z1
b_Z1 = b0_Z1 + b1_Z1*incidence
sig_Z1 = a_Z1*speed**(b_Z1)

# Z2
a0_Z2 = Z2_p[0]
a1_Z2 = Z2_p[1]
a2_Z2 = Z2_p[2]
b0_Z2 = Z2_p[3]
b1_Z2 = Z2_p[4]
b2_Z2 = Z2_p[5]

a_Z2 = a0_Z2 + a1_Z2*incidence + a2_Z2*incidence**2
b_Z2 = b0_Z2 + b1_Z2*incidence + b2_Z2*incidence**2
sig_Z2 = a_Z2*speed**(b_Z2)

c0, c1 = Final_p[0], Final_p[1]
c2, c3 = Final_p[2], Final_p[3]
sigmoid_Z1 = 1 / (1 + np.exp(-c0*(speed-c1)))
sigmoid_Z2 = 1 / (1 + np.exp(-c2*(speed-c3)))
sig_Final = sig_Z1 * sigmoid_Z1 + sig_Z2 * sigmoid_Z2
return sig_Final


@GmfModel.register(wspd_range=[3., 80.], pol='VH', units='linear')
def gmf_rcm_noaa(incidence, speed, phi=None):
"""
Radarsat Consteallation Mission VH GMF : relation between sigma0, incidence and windspeed.
From NOAA, named gmf_rcm_v1.
Parameters
----------
incidence: xarray.DataArray
incidence angle [deg]
speed: xarray.DataArray
wind speed [m/s]
Returns
-------
sigma0: xarray.DataArray
linear sigma0
"""
# constants params

Z1_p = np.array(
[2.2309436836414871e-12, 8.3374911282878728, -0.033443488982800210])
Z2_p = np.array([7.7945050373193260e-05, -2.4425748662769216e-06, 2.7625550632547159e-08,
1.2524896108831316, 0.019203092214131894, -0.00028408046502692580])
Final_p = np.array([-0.34498737004629487, 12.558975188752012,
0.12713502524515713, 4.2806865431046752])

# Z1

a0_Z1 = Z1_p[0]
b0_Z1 = Z1_p[1]
b1_Z1 = Z1_p[2]
a_Z1 = a0_Z1
b_Z1 = b0_Z1 + b1_Z1*incidence
sig_Z1 = a_Z1*speed**(b_Z1)
#Z2

a_Z1 = a0_Z1
b_Z1 = b0_Z1 + b1_Z1*incidence
sig_Z1 = a_Z1*speed**(b_Z1)

# Z2
a0_Z2 = Z2_p[0]
a1_Z2 = Z2_p[1]
a2_Z2 = Z2_p[2]
b0_Z2 = Z2_p[3]
b1_Z2 = Z2_p[4]
b2_Z2 = Z2_p[5]
a_Z2 = a0_Z2 + a1_Z2*incidence + a2_Z2*incidence**2

a_Z2 = a0_Z2 + a1_Z2*incidence + a2_Z2*incidence**2
b_Z2 = b0_Z2 + b1_Z2*incidence + b2_Z2*incidence**2
sig_Z2 = a_Z2*speed**(b_Z2)

c0,c1 = Final_p[0],Final_p[1]
c2,c3 = Final_p[2],Final_p[3]
c0, c1 = Final_p[0], Final_p[1]
c2, c3 = Final_p[2], Final_p[3]
sigmoid_Z1 = 1 / (1 + np.exp(-c0*(speed-c1)))
sigmoid_Z2 = 1 / (1 + np.exp(-c2*(speed-c3)))
sig_Final = sig_Z1 * sigmoid_Z1 + sig_Z2 * sigmoid_Z2
sig_Final = sig_Z1 * sigmoid_Z1 + sig_Z2 * sigmoid_Z2
return sig_Final
Loading

0 comments on commit 7b685bb

Please sign in to comment.