# # 大气延迟校正工具计算 # 参考：https://github.com/insarlab/PyAPS # from logHandler import log_info,log_err import logging import os import sys import scipy import numpy as np import scipy.interpolate as si from scipy import interpolate from interpolate import bilinearInterplor import netCDF4 as Nc import scipy.interpolate as intp import scipy.integrate as intg from ImageHandle import ImageHandler,ReprojectImages2,scatter2grid from osgeo import gdal,gdalconst from tqdm import tqdm from Process import PreProcess as pp from matplotlib import pyplot as plt logger = logging.getLogger("mylog") # 声明常数 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 class Gacosdelay: def __init__(self,output) -> None: self.wvl=None self.inc=None self.dlos=None self.lat=None self.lon=None self.dout=None self.pha=None self.output=output # def load_npy(self,output): # # 加载数据 # log_info("="*50) # log_info("加载预处理结果") # log_info("{}:{}".format(os.path.join(output,'lon.npy'),os.path.join(output,'lon.npy'))) # log_info("{}:{}".format(os.path.join(output,'lat.npy'),os.path.join(output,'lat.npy'))) # log_info("{}:{}".format(os.path.join(output,'wvl.npy'),os.path.join(output,'wvl.npy'))) # log_info("{}:{}".format(os.path.join(output,'inc.npy'),os.path.join(output,'inc.npy'))) # log_info("{}:{}".format(os.path.join(output,'dlos.tiff'),os.path.join(output,'dlos.tiff'))) # log_info("{}:{}".format(os.path.join(output,'dlosout.tiff'),os.path.join(output,'dlosout.tiff'))) # log_info("{}:{}".format(os.path.join(output,'pha.tiff'),os.path.join(output,'pha.tiff'))) # log_info("="*50) # self.lon=np.load(os.path.join(output,'lon.npy')) # self.lat=np.load(os.path.join(output,'lat.npy')) # self.wvl=np.load(os.path.join(output,'wvl.npy')) # self.inc=np.load(os.path.join(output,'inc.npy')) # self.dlos=ImageHandler.get_band_array(os.path.join(output,'dlos.tiff'),1) # self.dout=ImageHandler.get_band_array(os.path.join(output,'dlosout.tiff'),1) # self.pha=ImageHandler.get_band_array(os.path.join(output,'pha.tiff'),1) def save(self): # 保存成临时结果 np.save(os.path.join(self.output,'lon.npy'),np.array(self.lon)) np.save(os.path.join(self.output,'lat.npy'),np.array(self.lat)) np.save(os.path.join(self.output,'wvl.npy'),np.array(self.wvl)) np.save(os.path.join(self.output,'inc.npy'),np.array(self.inc)) ImageHandler.write_imgArray(os.path.join(self.output,'dlos.tiff'),self.dlos) ImageHandler.write_imgArray(os.path.join(self.output,'dlosout.tiff'),self.dout) ImageHandler.write_imgArray(os.path.join(self.output,'pha.tiff'),self.pha) 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 set_lon(self,reflect_lon_path): self.lon=ImageHandler.get_band_array(reflect_lon_path,1) def set_lat(self,reflect_lat_path): self.lat=ImageHandler.get_band_array(reflect_lat_path,1) def set_gacos_ztd(self,gacos_ztd_path,ztd_range_array_path): logger.info("对gacos ztd进行插值") im_geotrans = ImageHandler.get_dataset(gacos_ztd_path).GetGeoTransform() # 仿射矩阵 inv_gt=gdal.InvGeoTransform(im_geotrans) sample_in_tiff_x=inv_gt[0]+inv_gt[1]*self.lon+inv_gt[2]*self.lat # x sample_in_tiff_y=inv_gt[3]+inv_gt[4]*self.lon+inv_gt[5]*self.lat # y shape_range=sample_in_tiff_x.shape sample_in_tiff_x=sample_in_tiff_x.reshape(-1) sample_in_tiff_y=sample_in_tiff_y.reshape(-1) ztd_array = ImageHandler.get_band_array(gacos_ztd_path, 1) # 获取dem ztd_range_array=bilinearInterplor(ztd_array.astype(np.float64),sample_in_tiff_x.astype(np.float64),sample_in_tiff_y.astype(np.float64)) ztd_range_array=ztd_range_array.reshape(shape_range) ImageHandler.write_imgArray(ztd_range_array_path,ztd_range_array) self.dlos=ztd_range_array def get_delay(self): pha=self.dlos*np.pi*4.0/(np.cos(self.inc*np.pi/180.0)*self.wvl) self.dout=pha self.pha= -1*pha return self.pha def geo2range_dem(dem_path,reflect_lon_path,reflect_lat_path,out_range_dem_path ): lon_table=ImageHandler.get_band_array(reflect_lon_path,1) lat_table=ImageHandler.get_band_array(reflect_lat_path,1) im_geotrans = ImageHandler.get_dataset(dem_path).GetGeoTransform() # 仿射矩阵 inv_gt=gdal.InvGeoTransform(im_geotrans) sample_in_tiff_x=inv_gt[0]+inv_gt[1]*lon_table+inv_gt[2]*lat_table # x sample_in_tiff_y=inv_gt[3]+inv_gt[4]*lon_table+inv_gt[5]*lat_table # y shape_range=sample_in_tiff_x.shape sample_in_tiff_x=sample_in_tiff_x.reshape(-1) sample_in_tiff_y=sample_in_tiff_y.reshape(-1) dem_array = ImageHandler.get_band_array(dem_path, 1) # 获取dem dem_array=bilinearInterplor(dem_array.astype(np.float64),sample_in_tiff_x.astype(np.float64),sample_in_tiff_y.astype(np.float64)) dem_array=dem_array.reshape(shape_range) ImageHandler.write_imgArray(out_range_dem_path,dem_array) return True 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 =data - tropo shp=data.shape data=data-data[int(shp[0]/2),int(shp[1]/2)] # 校正 ImageHandler.write_imgArray(cor_dis_file,data) plt.matshow(data) plt.colorbar() plt.savefig(cor_dis_file.replace(".tiff",'.jpg').replace(".tif",'.jpg')) 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) plt.matshow(pha) plt.colorbar() plt.savefig(out_pha_path.replace(".tiff",'.jpg').replace(".tif",'.jpg')) 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: 从气象数据中读取数组 """ logger.info("="*50) logger.info("读取EAR(.nc) 文件：{}".format(nc_file)) logger.info("范围参数S,N,W,E:{}".format(str(Extent))) 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 logger.info("level:\n{}".format(str(self.level))) logger.info("lat_array:\n{}".format(str(self.lat_array))) logger.info("lon_array:\n{}".format(str(self.lon_array))) logger.info("r:\n{}".format(str(self.r))) logger.info("t:\n{}".format(str(self.t))) logger.info("z:\n{}".format(str(self.z))) # (1,level,lat,lon) except Exception as e : log_err(e) logger.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): logger.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() logger.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): logger.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 ## Save into files for plotting (test for era5 interpolation) #out_dir = '/home/angel/Tests/test_era5interpolation/Outputs' #suffix = '20141023_20141210_era5' #np.savetxt(f'{out_dir}/alt_{suffix}.txt', gph, fmt=' '.join(['%s']*380)) # altitude (of lvls) #np.savetxt(f'{out_dir}/lvls_{suffix}.txt', lvls.T, fmt='%s') # pressure #np.savetxt(f'{out_dir}/tmp_{suffix}.txt', tmp, fmt=' '.join(['%s']*380)) # temperature #np.savetxt(f'{out_dir}/vpr_{suffix}.txt', vpr, fmt=' '.join(['%s']*380)) # vapor #np.savetxt(f'{out_dir}/Alti_{suffix}.txt', hgt.T, fmt='%s') # interpo altitude #np.savetxt(f'{out_dir}/Presi_{suffix}.txt', Presi.T, fmt=' '.join(['%s']*380)) # interpo pressure #np.savetxt(f'{out_dir}/Tempi_{suffix}.txt', Tempi.T, fmt=' '.join(['%s']*380)) # interpo temperature #np.savetxt(f'{out_dir}/Vpri_{suffix}.txt', Vpri.T, fmt=' '.join(['%s']*380)) # interpo vapor 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 class tropo_atom_err: def __init__(self) -> None: pass def get_range_ztd_delay(self,ztddelay,reflect_table_path,geo2range=True): pass def get_delay(self,ztddelay,inc_path): pass def ztd_fromGeo2Range(self,in_ztd_file,reflect_table_path,out_range_ztd_path,work_space,geoflectTable=False,ifg_path=""): """ 将ztd 由地距空间映射到斜距空间中 Args: ztd_file (_type_): _description_ reflect_table_path (_type_): _description_ # lon，lat；rangeCoord_table，azimuthCoord_table out_range_ztd_path (_type_): _description_ geoflectTable (bool, optional): _description_. Defaults to False. """ ztd_array=None if geoflectTable: # 重采样到临时空间 temp_out_ztd_path=os.path.join(work_space,"ztd_geo_temp.tiff") if os.path.exists(temp_out_ztd_path): os.remove(temp_out_ztd_path) ReprojectImages2(in_ztd_file,reflect_table_path,temp_out_ztd_path) ztd_array = ImageHandler.get_band_array(temp_out_ztd_path, 1).reshape(-1) # 获取大气延迟表ztd # 根据映射表重新采样并生成坐标 rangeCoord_table=ImageHandler.get_band_array(reflect_table_path,1).reshape(-1) azimuthCoord_table=ImageHandler.get_band_array(reflect_table_path,2).reshape(-1) nncols = ImageHandler().get_img_width(ifg_path) nnrows = ImageHandler().get_img_height(ifg_path) mask=(rangeCoord_table0)*(azimuthCoord_table>0)*(azimuthCoord_table0 x_mask=np.where(x_mask==1)[0] x_arr=x_arr[x_mask] y_arr=y_arr[x_mask] lon_arr=lon_arr[x_mask] lat_arr=lat_arr[x_mask] dem_arr=dem_arr[x_mask] logger.info("插值得到 斜距->地距 lon") grid_lon_arr=scatter2grid(x_arr,y_arr,lon_arr,grid_x,grid_y) logger.info("插值得到 斜距->地距 lat") grid_lat_arr=scatter2grid(x_arr,y_arr,lat_arr,grid_x,grid_y) logger.info("插值得到 斜距->地距 dem") grid_dem_arr=scatter2grid(x_arr,y_arr,dem_arr,grid_x,grid_y) logger.info("插值得到 斜距->地距 结束") ImageHandler.write_imgArray(out_range_x_tiff_path,grid_lon_arr) ImageHandler.write_imgArray(out_range_y_tiff_path,grid_lat_arr) ImageHandler.write_imgArray(out_range_dem_path,grid_dem_arr) return True def ztd_fromGeo2Range(self,in_ztd_file,reflect_table_path,out_range_ztd_path,work_space,geoflectTable=False,ifg_path=""): """ 将ztd 由地距空间映射到斜距空间中 Args: ztd_file (_type_): _description_ reflect_table_path (_type_): _description_ # lon，lat；rangeCoord_table，azimuthCoord_table out_range_ztd_path (_type_): _description_ geoflectTable (bool, optional): _description_. Defaults to False. """ if geoflectTable: # 重采样到临时空间 temp_out_ztd_path=os.path.join(work_space,"ztd_geo_temp.tiff") if os.path.exists(temp_out_ztd_path): os.remove(temp_out_ztd_path) ReprojectImages2(in_ztd_file,reflect_table_path,temp_out_ztd_path) ztd_array = ImageHandler.get_band_array(temp_out_ztd_path, 1).reshape(-1) # 获取大气延迟表ztd # 根据映射表重新采样并生成坐标 rangeCoord_table=ImageHandler.get_band_array(reflect_table_path,1).reshape(-1) azimuthCoord_table=ImageHandler.get_band_array(reflect_table_path,2).reshape(-1) nncols = ImageHandler().get_img_width(ifg_path) nnrows = ImageHandler().get_img_height(ifg_path) mask=(rangeCoord_table0)*(azimuthCoord_table>0)*(azimuthCoord_table