499 lines
19 KiB
Python
499 lines
19 KiB
Python
#
|
||
# 大气延迟校正工具计算
|
||
# 参考：https://github.com/insarlab/PyAPS
|
||
#
|
||
|
||
import os
|
||
import numpy as np
|
||
import scipy.interpolate as si
|
||
import netCDF4 as Nc
|
||
import scipy.interpolate as intp
|
||
import scipy.integrate as intg
|
||
from tool.algorithm.image.ImageHandle import ImageHandler
|
||
from tqdm import tqdm
|
||
|
||
# 声明常数
|
||
constdict = {}
|
||
k1 = 0.776 #(K/Pa)
|
||
k2 = 0.716 #(K/Pa)
|
||
k3 = 3750 #(K^2.Pa)
|
||
g = 9.81 #(m/s^2)
|
||
Rd = 287.05 #(J/Kg/K)
|
||
Rv = 461.495 #(J/Kg/K)
|
||
mma = 29.97 #(g/mol)
|
||
mmH = 2.0158 #(g/mol)
|
||
mmO = 16.0 #(g/mol)
|
||
Rho = 1000.0 #(kg/m^3)
|
||
|
||
a1w = 611.21 # hPa
|
||
a3w = 17.502 #
|
||
a4w = 32.19 # K
|
||
a1i = 611.21 # hPa
|
||
a3i = 22.587 #
|
||
a4i = -0.7 # K
|
||
T3 = 273.16 # K
|
||
Ti = 250.16 # K
|
||
nhgt = 300 # Number of levels for interpolation (OLD 151)
|
||
minAlt = -200.0
|
||
maxAlt = 50000.0
|
||
minAltP = -200.0
|
||
|
||
|
||
def correct_single_ifgram(ifg_path, phas_path, cor_dis_file):
|
||
data=ImageHandler.get_band_array(ifg_path,1)
|
||
tropo=ImageHandler.get_band_array(phas_path,1)
|
||
data =tropo - np.angle(data)
|
||
return data
|
||
|
||
def createpha(master_pha,aux_pha,wavelen,out_pha_path):
|
||
pha=master_pha-aux_pha
|
||
#pha=pha*-4. * np.pi / float( wavelen)
|
||
ImageHandler.write_imgArray(out_pha_path,pha)
|
||
return pha
|
||
|
||
class ERA:
|
||
""" 用于存储从ear 解析得到的参数"""
|
||
def __init__(self) -> None:
|
||
# 定义计算变量
|
||
self.lon_array=None
|
||
self.lat_array=None
|
||
self.level=None
|
||
self.t=None
|
||
self.r=None
|
||
self.z=None
|
||
pass
|
||
|
||
def read_nc(self,nc_file,Extent=None,level_list=None):
|
||
"""
|
||
function: 从气象数据中读取数组
|
||
|
||
"""
|
||
try:
|
||
nc_data_obj = Nc.Dataset(nc_file)
|
||
inform = nc_data_obj.variables
|
||
lon_array = inform['longitude'][:] # lon数组
|
||
lat_array = inform['latitude'][:] # lat数组
|
||
n_lat = lat_array.shape[0] # 经线列数量
|
||
n_lon = lon_array.shape[0] # 纬线行数量
|
||
level=np.array(inform["level"]) # pressure level in m 压力面
|
||
t = np.array(inform["t"]) # temp--k (1,37,latitude,longitude)
|
||
r = np.array(inform["r"]) # relative_humidity--% (1,37,latitude,longitude)
|
||
z = np.array(inform["z"]) # Geopotential--m**2 s**-2 (1,37,latitude,longitude)
|
||
|
||
if not Extent is None:
|
||
min_lat,max_lat,min_lon,max_lon=Extent
|
||
mask_lon= (lon_array>=min_lon)&(lon_array<=max_lon)
|
||
mask_lat= (lat_array>=min_lat)&(lat_array<=max_lat)
|
||
mask_lon=np.where(mask_lon==1)[0]
|
||
mask_lat=np.where(mask_lat==1)[0]
|
||
lat_array=lat_array[mask_lat]
|
||
lon_array=lat_array[mask_lon]
|
||
t=t[:,:,mask_lat,:][:,:,:,mask_lon]
|
||
r=r[:,:,mask_lat,:][:,:,:,mask_lon]
|
||
z=z[:,:,mask_lat,:][:,:,:,mask_lon]
|
||
if not level_list is None:
|
||
mask_level=[]
|
||
for l_item in level_list:
|
||
mask_level.append(np.where(level==l_item)[0][0])
|
||
level=level[mask_level]
|
||
t=t[:,mask_level,:,:]
|
||
r=r[:,mask_level,:,:]
|
||
z=z[:,mask_level,:,:]
|
||
# 根据
|
||
self.level=level
|
||
self.lat_array=lat_array
|
||
self.lon_array=lon_array
|
||
self.r=r
|
||
self.t=t
|
||
self.z=z
|
||
# log_info("level:\n{}".format(str(self.level)))
|
||
# log_info("lat_array:\n{}".format(str(self.lat_array)))
|
||
# log_info("lon_array:\n{}".format(str(self.lon_array)))
|
||
# log_info("r:\n{}".format(str(self.r)))
|
||
# log_info("t:\n{}".format(str(self.t)))
|
||
# log_info("z:\n{}".format(str(self.z))) # (1,level,lat,lon)
|
||
except Exception as e :
|
||
print(e)
|
||
# log_err(e)
|
||
# log_info("="*50)
|
||
return True
|
||
|
||
|
||
class AtomZtdDelay:
|
||
def __init__(self,atomdelay_path) -> None:
|
||
self.init_var()
|
||
self.const_var()
|
||
self.output=atomdelay_path
|
||
pass
|
||
def const_var(self):
|
||
self.lvls = np.array([1, 2, 3, 5, 7, 10, 20, 30, 50, 70, 100, 125,
|
||
150, 175,200, 225, 250, 300, 350, 400, 450,
|
||
500, 550, 600, 650, 700, 750, 775, 800, 825,
|
||
850, 875, 900, 925, 950, 975, 1000])
|
||
self.nlvls = len(self.lvls)
|
||
self.hgtscale = ((maxAlt - minAlt) / nhgt) / 0.703
|
||
self.bufspc = 1.2
|
||
self.alpha = Rv / Rd
|
||
|
||
def set_waveLen(self,wvl):
|
||
self.wvl=wvl
|
||
|
||
def set_incidenceAngle(self,inc_path):
|
||
local_angle_array = ImageHandler.get_band_array(inc_path, 1) # 获取局部入射角
|
||
self.inc=local_angle_array
|
||
|
||
def init_var(self):
|
||
self.deley_all=None
|
||
self.Dry_delay=None
|
||
self.Wet_delay=None
|
||
self.lon_arr=None
|
||
self.lat_arr=None
|
||
self.Pi=None
|
||
self.Ti=None
|
||
self.Vi=None
|
||
self.hgt=None
|
||
self.lvls=None
|
||
self.gph=None
|
||
self.tmp=None
|
||
self.vpr=None
|
||
self.dout=None
|
||
|
||
def load_npy(self,output):
|
||
# 加载数据
|
||
# log_info("="*50)
|
||
# log_info("加载预处理结果")
|
||
# log_info("{}:{}".format(os.path.join(output,'deley_all.npy'),os.path.exists(os.path.join(output,'deley_all.npy'))))
|
||
# log_info("{}:{}".format(os.path.join(output,'Dry_delay.npy'),os.path.exists(os.path.join(output,'Dry_delay.npy'))))
|
||
# log_info("{}:{}".format(os.path.join(output,'Wet_delay.npy'),os.path.exists(os.path.join(output,'Wet_delay.npy'))))
|
||
# log_info("{}:{}".format(os.path.join(output,'lon_arr.npy'),os.path.exists(os.path.join(output,'lon_arr.npy'))))
|
||
# log_info("{}:{}".format(os.path.join(output,'lat_arr.npy'),os.path.exists(os.path.join(output,'lat_arr.npy'))))
|
||
# log_info("{}:{}".format(os.path.join(output,'Pi.npy'),os.path.exists(os.path.join(output,'Pi.npy'))))
|
||
# log_info("{}:{}".format(os.path.join(output,'Ti.npy'),os.path.exists(os.path.join(output,'Ti.npy'))))
|
||
# log_info("{}:{}".format(os.path.join(output,'Vi.npy'),os.path.exists(os.path.join(output,'Vi.npy'))))
|
||
# log_info("{}:{}".format(os.path.join(output,'hgt.npy'),os.path.exists(os.path.join(output,'hgt.npy'))))
|
||
# log_info("{}:{}".format(os.path.join(output,'lvls.npy'),os.path.exists(os.path.join(output,'lvls.npy'))))
|
||
# log_info("{}:{}".format(os.path.join(output,'gph.npy'),os.path.exists(os.path.join(output,'gph.npy'))))
|
||
# log_info("{}:{}".format(os.path.join(output,'tmp.npy'),os.path.exists(os.path.join(output,'tmp.npy'))))
|
||
# log_info("{}:{}".format(os.path.join(output,'vpr.npy'),os.path.exists(os.path.join(output,'vpr.npy'))))
|
||
# log_info("{}:{}".format(os.path.join(output,'dout.npy'),os.path.exists(os.path.join(output,'dout.npy'))))
|
||
# log_info("="*50)
|
||
self.deley_all=np.load(os.path.join(output,'deley_all.npy'))
|
||
self.Dry_delay=np.load(os.path.join(output,'Dry_delay.npy'))
|
||
self.Wet_delay=np.load(os.path.join(output,'Wet_delay.npy'))
|
||
self.lon_arr=np.load(os.path.join(output,'lon_arr.npy'))
|
||
self.lat_arr=np.load(os.path.join(output,'lat_arr.npy'))
|
||
self.Pi=np.load(os.path.join(output,'Pi.npy'))
|
||
self.Ti=np.load(os.path.join(output,'Ti.npy'))
|
||
self.Vi=np.load(os.path.join(output,'Vi.npy'))
|
||
self.hgt=np.load(os.path.join(output,'hgt.npy'))
|
||
self.lvls=np.load(os.path.join(output,'lvls.npy')).tolist()
|
||
self.gph=np.load(os.path.join(output,'gph.npy'))
|
||
self.tmp=np.load(os.path.join(output,'tmp.npy'))
|
||
self.vpr=np.load(os.path.join(output,'vpr.npy'))
|
||
self.dout=np.load(os.path.join(output,'dlosout.npy'))
|
||
self.pha=np.load(os.path.join(output,'pha.npy'))
|
||
|
||
def save(self):
|
||
# 保存成临时结果
|
||
ImageHandler.write_imgArray(os.path.join(self.output,'deley_all.tiff'),self.deley_all)
|
||
# np.save(os.path.join(self.output,'Dry_delay.npy'),self.Dry_delay)
|
||
# np.save(os.path.join(self.output,'Wet_delay.npy'),self.Wet_delay)
|
||
# np.save(os.path.join(self.output,'lon_arr.npy'),np.array(self.lon_arr))
|
||
# np.save(os.path.join(self.output,'lat_arr.npy'),np.array(self.lat_arr))
|
||
# np.save(os.path.join(self.output,'Pi.npy'),self.Pi)
|
||
# np.save(os.path.join(self.output,'Ti.npy'),self.Ti)
|
||
# np.save(os.path.join(self.output,'Vi.npy'),self.Vi)
|
||
# np.save(os.path.join(self.output,'hgt.npy'),self.hgt)
|
||
# np.save(os.path.join(self.output,'lvls.npy'),np.array(self.lvls))
|
||
# np.save(os.path.join(self.output,'gph.npy'),self.gph)
|
||
# np.save(os.path.join(self.output,'tmp.npy'),self.tmp)
|
||
# np.save(os.path.join(self.output,'vpr.npy'),self.vpr)
|
||
ImageHandler.write_imgArray(os.path.join(self.output,'dlosout.tiff'),self.dout) # ztd
|
||
ImageHandler.write_imgArray(os.path.join(self.output,'pha.tiff'),self.pha)
|
||
|
||
def get_dem(self,dem_path):
|
||
dem_array=ImageHandler.get_band_array(dem_path,1)
|
||
self.dem=dem_array
|
||
return True
|
||
|
||
def get_lon_lat(self,reflect_lon_path,reflect_lat_path):
|
||
lon_table=ImageHandler.get_band_array(reflect_lon_path,1)
|
||
lat_table=ImageHandler.get_band_array(reflect_lat_path,1)
|
||
|
||
self.lon=lon_table
|
||
self.lat=lat_table
|
||
return True
|
||
|
||
def get_delay(self):
|
||
# log_info("计算相位")
|
||
minAltp=-200.0
|
||
|
||
lon_arr=self.lon
|
||
lat_arr=self.lat
|
||
# 根据行列号
|
||
min_lon=lon_arr.min()-self.bufspc
|
||
max_lon=lon_arr.max()+self.bufspc
|
||
min_lat=lat_arr.min()-self.bufspc
|
||
max_lat=lat_arr.max()+self.bufspc
|
||
# WGS [-180,-90,180,90]
|
||
# dem
|
||
shp=lon_arr.shape
|
||
lon_arr=lon_arr.reshape(-1)
|
||
if min_lon<0 and max_lon>0:
|
||
lon_arr[np.where(lon_arr<0)]=lon_arr[np.where(lon_arr<0)]+360
|
||
lon_arr=lon_arr.reshape(shp)
|
||
self.ny,self.nx=self.dem.shape
|
||
|
||
dout = np.zeros((self.ny, self.nx), dtype=np.float32)
|
||
|
||
intp_1d = si.interp1d(self.hgt, self.deley_all, kind='cubic', axis=-1)
|
||
|
||
self.dem[np.isnan(self.dem)] = minAltp
|
||
self.dem[self.dem < minAltp] = minAltp
|
||
self.mask=np.ones(self.dem.shape)
|
||
minH = np.max([np.nanmin(self.dem), self.hgt.min()])
|
||
maxH = int(np.nanmax(self.dem)) + 100.
|
||
kh = np.arange(minH,maxH)
|
||
self.Delfn_1m = intp_1d(kh)
|
||
self.alti = kh
|
||
# no reshape
|
||
|
||
Lonu = self.lon_arr[0,:]
|
||
Latu = self.lat_arr[:,0]
|
||
linearint = si.RegularGridInterpolator(
|
||
(Latu[::-1], Lonu,kh),
|
||
self.Delfn_1m[::-1,:,:],
|
||
method='linear',
|
||
bounds_error=False,
|
||
fill_value=0.0,
|
||
)
|
||
wvl=self.wvl
|
||
# Loop on the lines
|
||
pha=dout.copy()
|
||
# log_info("计算延迟值")
|
||
for m in tqdm(range(self.ny)):
|
||
|
||
# Get latitude and longitude arrays
|
||
lati = self.lat[m,:]*self.mask[m,:]
|
||
loni = self.lon[m,:]*self.mask[m,:]
|
||
|
||
# Remove negative values
|
||
|
||
# Remove NaN values
|
||
ii = np.where(np.isnan(lati))
|
||
jj = np.where(np.isnan(loni))
|
||
xx = np.union1d(ii,jj)
|
||
lati[xx]=0.0
|
||
loni[xx]=0.0
|
||
|
||
cinc = np.cos(self.inc[m,:]*np.pi/180.)
|
||
|
||
# Make the bilinear interpolation
|
||
D = self.dem[m,:]
|
||
val = linearint(np.vstack((lati, loni, D)).T)
|
||
val[xx] = np.nan
|
||
|
||
# save output
|
||
dout[m,:] = val
|
||
pha[m,:]=val*np.pi*4.0/(cinc*wvl)
|
||
|
||
self.dout=dout
|
||
self.pha= -1*pha
|
||
return pha
|
||
pass
|
||
|
||
def get_deley_dlf(self,nc_file):
|
||
# log_info("计算延迟路径")
|
||
[lon_array,lat_array, lvls, gph,tmp,vpr]=self.get_ear5_nc(nc_file)
|
||
lon_array,lat_array=np.meshgrid(lon_array,lat_array)
|
||
hgt = np.linspace(minAltP, gph.max().round(), nhgt)
|
||
[Pi,Ti,Vi]=self.intP2H(lvls,hgt,gph,tmp,vpr)
|
||
[Dry_delay,Wet_delay] = self.PTV2del(Pi,Ti,Vi,hgt)
|
||
deley_all = Dry_delay+Wet_delay
|
||
# 保存结果
|
||
self.deley_all=deley_all
|
||
self.Dry_delay=Dry_delay
|
||
self.Wet_delay=Wet_delay
|
||
self.lon_arr=lon_array
|
||
self.lat_arr=lat_array
|
||
self.Pi=Pi
|
||
self.Ti=Ti
|
||
self.Vi=Vi
|
||
self.hgt=hgt
|
||
self.lvls=lvls
|
||
self.gph=gph # 高度
|
||
self.tmp=tmp # 温度
|
||
self.vpr=vpr # 湿度
|
||
return True
|
||
|
||
def get_ear5_nc(self,nc_file,Extent=None):
|
||
era= ERA()
|
||
era.read_nc(nc_file,Extent,self.lvls)
|
||
lvls=self.lvls*100 # 转为 绝对压力
|
||
nlvls=len(era.level)
|
||
nlat=len(era.lat_array)
|
||
nlon=len(era.lon_array)
|
||
# 逐层分析
|
||
gph=1.0*era.z[0,:,:,:]/g
|
||
tmp=era.t[0,:,:,:]
|
||
vpr=era.r[0,:,:,:]
|
||
for i in range(nlvls):
|
||
esat = self.cc_era(tmp[i,:,:])
|
||
vpr[i,:,:]=(vpr[i,:,:]/100.0)*esat
|
||
#
|
||
return era.lon_array,era.lat_array, lvls, gph,tmp,vpr
|
||
|
||
def intP2H(self,lvls,hgt,gph,tmp,vpr):
|
||
# Interpolate pressure, temperature and Humidity of hgt
|
||
#minAlt = minAlt #Hardcoded parameter. -200.0
|
||
maxAlt_2 = gph.max().round()
|
||
|
||
nlat = gph.shape[1] #Number of stations
|
||
nlon = gph.shape[2]
|
||
nhgt = len(hgt) #Number of height points
|
||
Presi = np.zeros((nlat,nlon,nhgt))
|
||
Tempi = np.zeros((nlat,nlon,nhgt))
|
||
Vpri = np.zeros((nlat,nlon,nhgt))
|
||
|
||
for i in range(nlat):
|
||
for j in range(nlon):
|
||
temp = gph[:,i,j] #Obtaining height values
|
||
hx = temp.copy()
|
||
sFlag = False
|
||
eFlag = False
|
||
#Add point at start
|
||
if (hx.min() > minAlt): # -200.0
|
||
sFlag = True
|
||
#changed from 1 to 0 (-1 should also work), CL
|
||
hx = np.concatenate((hx,[minAlt-1]),axis = 0)
|
||
|
||
#Add point at end
|
||
if (hx.max() < maxAlt_2):
|
||
eFlag = True
|
||
#changed from 1 to 0 (-1 should also work), CL
|
||
hx = np.concatenate(([maxAlt_2+1],hx),axis=0)
|
||
|
||
#Splines needs monotonically increasing.
|
||
hx = -hx
|
||
|
||
#Interpolating pressure values
|
||
hy = lvls.copy()
|
||
if (sFlag == True):
|
||
val = hy[-1] +(hx[-1] - hx[-2])* (hy[-1] - hy[-2])/(hx[-2]-hx[-3])
|
||
#changed from 1 to 0 (-1 should also work), CL
|
||
hy = np.concatenate((hy,[val]),axis=0)
|
||
|
||
if (eFlag == True):
|
||
val = hy[0] - (hx[0] - hx[1]) * (hy[0] - hy[1])/(hx[1]-hx[2])
|
||
#changed from 1 to 0 (-1 should also work), CL
|
||
hy = np.concatenate(([val],hy),axis=0)
|
||
|
||
tck = intp.interp1d(hx,hy,kind='cubic')
|
||
#Again negative for consistency with hx
|
||
temp = tck(-hgt)
|
||
Presi[i,j,:] = temp.copy()
|
||
del temp
|
||
|
||
#Interpolating temperature
|
||
temp = tmp[:,i,j]
|
||
hy = temp.copy()
|
||
if (sFlag == True):
|
||
val = hy[-1] +(hx[-1] - hx[-2])* (hy[-1] - hy[-2])/(hx[-2]-hx[-3])
|
||
#changed from 1 to 0 (-1 should also work), CL
|
||
hy = np.concatenate((hy,[val]),axis=0)
|
||
|
||
if (eFlag == True):
|
||
val = hy[0] - (hx[0] - hx[1]) * (hy[0] - hy[1])/(hx[1]-hx[2])
|
||
#changed from 1 to 0 (-1 should also work), CL
|
||
hy = np.concatenate(([val],hy),axis=0)
|
||
|
||
tck = intp.interp1d(hx,hy,kind='cubic')
|
||
temp = tck(-hgt)
|
||
Tempi[i,j,:] = temp.copy()
|
||
del temp
|
||
|
||
#Interpolating vapor pressure
|
||
temp = vpr[:,i,j]
|
||
hy = temp.copy()
|
||
if (sFlag == True):
|
||
val = hy[-1] +(hx[-1] - hx[-2])* (hy[-1] - hy[-2])/(hx[-2]-hx[-3])
|
||
#changed from 1 to 0 (-1 should also work), CL
|
||
hy = np.concatenate((hy,[val]),axis=0)
|
||
|
||
if (eFlag == True):
|
||
val = hy[0] - (hx[0] - hx[1]) * (hy[0] - hy[1])/(hx[1]-hx[2])
|
||
#changed from 1 to 0 (-1 should also work), CL
|
||
hy = np.concatenate(([val],hy),axis=0)
|
||
|
||
tck = intp.interp1d(hx,hy,kind='cubic')
|
||
temp = tck(-hgt)
|
||
Vpri[i,j,:] = temp.copy()
|
||
del temp
|
||
|
||
return Presi,Tempi,Vpri
|
||
|
||
def PTV2del(self,Presi,Tempi,Vpri,hgt):
|
||
'''Computes the delay function given Pressure, Temperature and Vapor pressure.
|
||
|
||
Args:
|
||
* Presi (np.ndarray) : Pressure at height levels.
|
||
* Tempi (np.ndarray) : Temperature at height levels.
|
||
* Vpri (np.ndarray) : Vapor pressure at height levels.
|
||
* hgt (np.ndarray) : Height levels.
|
||
* cdict (np.ndarray) : Dictionary of constants.
|
||
|
||
Returns:
|
||
* DDry2 (np.ndarray) : Dry component of atmospheric delay.
|
||
* DWet2 (np.ndarray) : Wet component of atmospheric delay.
|
||
|
||
.. note::
|
||
Computes refractive index at each altitude and integrates the delay using cumtrapz.'''
|
||
nhgt = len(hgt) #Number of height points
|
||
nlat = Presi.shape[0] #Number of stations
|
||
nlon = Presi.shape[1]
|
||
WonT = Vpri/Tempi
|
||
WonT2 = WonT/Tempi
|
||
|
||
|
||
#Dry delay
|
||
DDry2 = np.zeros((nlat,nlon,nhgt))
|
||
DDry2[:,:,:] = k1*Rd*(Presi[:,:,:] - Presi[:,:,-1][:,:,np.newaxis])*1.0e-6/g
|
||
|
||
#Wet delay
|
||
S1 = intg.cumtrapz(WonT,x=hgt,axis=-1)
|
||
val = 2*S1[:,:,-1]-S1[:,:,-2]
|
||
val = val[:,:,None]
|
||
S1 = np.concatenate((S1,val),axis=-1)
|
||
del WonT
|
||
|
||
S2 = intg.cumtrapz(WonT2,x=hgt,axis=-1)
|
||
val = 2*S2[:,:,-1]-S2[:,:,-2]
|
||
val = val[:,:,None]
|
||
S2 = np.concatenate((S2,val),axis=-1)
|
||
DWet2 = -1.0e-6*((k2-k1*Rd/Rv)*S1+k3*S2)
|
||
|
||
for i in range(nlat):
|
||
for j in range(nlon):
|
||
DWet2[i,j,:] = DWet2[i,j,:] - DWet2[i,j,-1]
|
||
|
||
return DDry2,DWet2
|
||
|
||
def cc_era(self,tmp):
|
||
shp=tmp.shape
|
||
tmp=tmp.reshape(-1)
|
||
esatw = a1w*np.exp(a3w*(tmp-T3)/(tmp-a4w))
|
||
esati = a1i*np.exp(a3i*(tmp-T3)/(tmp-a4i))
|
||
result = np.zeros(esati.shape)
|
||
mask_T3=tmp>=T3
|
||
mask_Ti=tmp<=Ti
|
||
mask_wgt=(mask_T3+mask_Ti)==0
|
||
mask_T3=np.where(mask_T3==1)
|
||
mask_Ti=np.where(mask_Ti==1)
|
||
mask_wgt=np.where(mask_wgt==1)
|
||
result[mask_T3]=esatw[mask_T3]
|
||
result[mask_Ti]=esati[mask_Ti]
|
||
result[mask_wgt]=esati[mask_wgt]+(esatw[mask_wgt]-esati[mask_wgt])*np.square((tmp[mask_wgt]-Ti)/(T3-Ti))
|
||
result=result.reshape(shp)
|
||
return result
|
||
|