129 lines
4.8 KiB
Cython
129 lines
4.8 KiB
Cython
# -*- coding: utf-8 -*-
|
|
"""
|
|
Created on Tue Jun 4 14:59:54 2013
|
|
|
|
@author: Sat Kumar Tomer
|
|
@email: satkumartomer@gmail.com
|
|
@website: www.ambhas.com
|
|
|
|
"""
|
|
cimport cython # 必须导入
|
|
import numpy as np##必须为c类型和python类型的数据都申明一个np
|
|
cimport numpy as np # 必须为c类型和python类型的数据都申明一个np
|
|
from libc.math cimport pi
|
|
from scipy.optimize import fmin
|
|
|
|
cpdef np.ndarray[double,ndim=1] inverse_oh2004(double sigma0vvdB,double sigma0hhdB,double sigma0hvdB,double theta,double f):
|
|
"""
|
|
sigma0vvdB = -14.1 dB
|
|
sigma0hhdB = -16.0
|
|
sigma0hvdB = -26.5
|
|
theta = 35. 角度值 ##Incidence angle
|
|
f = 5.0 ##GHz
|
|
"""
|
|
#print("--------------------------------------------------------\n")
|
|
cdef np.ndarray[double,ndim=1] result=np.ones((2))
|
|
result[0]=np.nan
|
|
result[1]=np.nan
|
|
#print("*************设置为nan****************")
|
|
#print(sigma0vvdB,sigma0hhdB,sigma0hvdB,theta,f)
|
|
cdef double k = 2*3.1415926*f/0.299792458; #calculate the wave number
|
|
cdef double theta_rad = theta*3.1415926/180; #represent angle in radians
|
|
|
|
cdef double sigma_0_vv = np.power(10.,(sigma0vvdB/10.)) #%represent data in linear scale
|
|
cdef double sigma_0_hh = np.power(10.,(sigma0hhdB/10.))
|
|
cdef double sigma_0_hv = np.power(10.,(sigma0hvdB/10.))
|
|
|
|
if sigma_0_vv==0:
|
|
#print("***********sigma_0_vv==0*************")
|
|
return result
|
|
cdef double p = sigma_0_hh / sigma_0_vv; #calculate the p-ratio
|
|
cdef double q = sigma_0_hv / sigma_0_vv; #calculate the q-ratio
|
|
|
|
|
|
cdef np.ndarray[double,ndim=1] mv0 = np.arange(0.05,0.9,0.01) # set Gamma0 range of values (fine increments)
|
|
|
|
## First estimates s1 and mv1
|
|
cdef np.ndarray[double,ndim=1] ks = ((-1.)*3.125*np.log(1 - sigma_0_hv/(0.11 * mv0**0.7 * (np.cos(theta_rad))**2.2)))**0.556
|
|
cdef np.ndarray[double,ndim=1] err = (1. - (2.*theta_rad/np.pi)**(0.35*mv0**(-0.65)) * np.exp(-0.4 * ks**1.4))-p
|
|
cdef np.ndarray[double,ndim=1] abs_err = np.abs(err);
|
|
cdef double min_err = np.nanmin(abs_err); #find the value of minimum error
|
|
#print(np.where(abs_err == min_err)[0].shape)
|
|
if min_err==np.nan or np.max(np.where(abs_err == min_err)[0].shape)==0 :
|
|
#print("***************min_err==np.nan or np.max(np.where(abs_err == min_err)[0].shape)==0")
|
|
return result
|
|
cdef double mv1 = mv0[np.where(abs_err == min_err)[0][0]]
|
|
cdef double temp_ks1=1. - sigma_0_hv/(0.11 * mv1**0.7 * (np.cos(theta_rad))**2.2)
|
|
if temp_ks1<0:
|
|
#print("*********************temp_ks1<0")
|
|
return result
|
|
cdef double ks1 = ((-1)*3.125*np.log(temp_ks1))**0.556
|
|
cdef double s1 = ks1/k
|
|
|
|
## Second estimate s2 and mv2
|
|
cdef double ks2 = (np.log(1-(q/(0.095 * (0.13 + np.sin(1.5*theta_rad))**1.4))) /(-1.3))**(10./9.)
|
|
cdef double s2 = ks2/k
|
|
cdef double mv2 =0.
|
|
cdef double yy =0.
|
|
cdef double xx = (1-p)/np.exp(-0.4 * ks2**1.4)
|
|
if xx<=0:
|
|
mv2 =0.
|
|
else:
|
|
yy = np.log(xx)/(0.35*np.log(2*theta_rad/np.pi))
|
|
mv2=np.power(yy,-100.0/65)
|
|
|
|
## Third estimate mv3
|
|
cdef double mv3 = ((sigma_0_hv/(1 - np.exp(-0.32 * ks2**1.8)))/(0.11 * np.cos(theta_rad)**2.2))**(1/0.7)
|
|
## weighted average s and mv-------------------------------------
|
|
#print("q:\t",q)
|
|
#print("k:\t",k)
|
|
#print("ks1:\t",ks1)
|
|
#print("ks2:\t",ks2)
|
|
#print("theta_rad:\t",theta_rad)
|
|
|
|
cdef double sf = (s1 + 0.25*s2)/(1+0.25)
|
|
cdef double mvf = (mv1+mv2+mv3)/3
|
|
|
|
result[0]=mvf*1.0
|
|
result[1]=sf*1.0
|
|
#print("mv1:\t",mv1)
|
|
#print("mv2:\t",mv2)
|
|
#print("mv3:\t",mv3)
|
|
#print("s1:\t",s1)
|
|
#print("s2:\t",s2)
|
|
#print("Estimated volumetric soil moisture: ",result[0])
|
|
#print("Estimated rms height s (m): ",result[1])
|
|
#print("\nend\n")
|
|
return result
|
|
|
|
cpdef double lamda2freq(double lamda):
|
|
return 299792458.0/lamda
|
|
|
|
cpdef double freq2lamda(double freq):
|
|
return 299792458.0/freq
|
|
|
|
# double sigma0vvdB,double sigma0hhdB,double sigma0hvdB,double theta,double f
|
|
cpdef int retrieve_oh2004_main(int n,np.ndarray[double,ndim=1] mv,np.ndarray[double,ndim=1] h,np.ndarray[int,ndim=1] mask,np.ndarray[double,ndim=1] sigma0vvdB,np.ndarray[double,ndim=1] sigma0hhdB,np.ndarray[double,ndim=1] sigma0hvdB, np.ndarray[double,ndim=1] vh, np.ndarray[double,ndim=1] theta,double f):
|
|
cdef int i=0;
|
|
cdef np.ndarray[double,ndim=1] result;
|
|
while i<n:
|
|
if mask[i]<0.5:
|
|
mv[i]=np.nan
|
|
h[i] =np.nan
|
|
else:
|
|
#print(i)
|
|
##print(sigma0vvdB[i], sigma0hhdB[i],sigma0hvdB[i], theta[i], f)
|
|
result= inverse_oh2004(sigma0vvdB[i], sigma0hhdB[i],sigma0hvdB[i], theta[i], f)
|
|
##print(result)
|
|
mv[i]=result[0]
|
|
h[i] =result[1]
|
|
##print(mv[i],h[i])
|
|
##print(result[0],result[1])
|
|
i=i+1
|
|
return 1
|
|
|
|
|
|
|
|
|
|
|