884 lines
37 KiB
Python
884 lines
37 KiB
Python
|
#
|
|||
|
|
# 大气延迟校正工具计算
|
|||
|
|
# 参考：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_table<nncols)*(rangeCoord_table>0)*(azimuthCoord_table>0)*(azimuthCoord_table<nnrows) # 注意0区域
|
|||
|
|
x=rangeCoord_table[np.where(mask==1)].reshape(-1)
|
|||
|
|
y=azimuthCoord_table[np.where(mask==1)].reshape(-1)
|
|||
|
|
z=ztd_array[np.where(mask==1)].reshape(-1)
|
|||
|
|
nx,ny=np.meshgrid(np.arange(nncols),np.arange(nnrows))
|
|||
|
|
nx=nx.reshape(-1)
|
|||
|
|
ny=ny.reshape(-1)
|
|||
|
|
ztd_range_array = scipy.interpolate.griddata((y,x),z,(ny,nx))
|
|||
|
|
ztd_range_array=ztd_range_array.reshape(nnrows,nncols)
|
|||
|
|
ImageHandler.write_imgArray(out_range_ztd_path,ztd_range_array)
|
|||
|
|
else:
|
|||
|
|
lon_table=ImageHandler.get_band_array(reflect_table_path,1)
|
|||
|
|
lat_table=ImageHandler.get_band_array(reflect_table_path,2)
|
|||
|
|
#
|
|||
|
|
im_geotrans = ImageHandler.get_dataset(ztd_array).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)
|
|||
|
|
# 斜距DLOS
|
|||
|
|
ztd_array = ImageHandler.get_band_array(in_ztd_file, 1) # 获取大气延迟表ztd
|
|||
|
|
logger.info("project ztd from geo to sar :{}".format(out_range_ztd_path))
|
|||
|
|
ztd_array=bilinearInterplor(ztd_array,sample_in_tiff_x,sample_in_tiff_y)
|
|||
|
|
ztd_array=ztd_array.reshap(shape_range)
|
|||
|
|
ImageHandler.write_imgArray(out_range_ztd_path,ztd_array)
|
|||
|
|
|
|||
|
|
|
|||
|
|
def geo2range(in_geo_x_tiff_path,in_geo_y_tiff_path,ref_ifg_path,dem_path,out_range_x_tiff_path,out_range_y_tiff_path,out_range_dem_path):
|
|||
|
|
geo_height=ImageHandler.get_img_height(in_geo_x_tiff_path)
|
|||
|
|
geo_width=ImageHandler.get_img_width(in_geo_x_tiff_path)
|
|||
|
|
range_height=ImageHandler.get_img_height(ref_ifg_path)
|
|||
|
|
range_width=ImageHandler.get_img_width(ref_ifg_path)
|
|||
|
|
|
|||
|
|
geotrans=ImageHandler.get_geotransform(in_geo_x_tiff_path)
|
|||
|
|
geo_x_arr,geo_y_arr=np.meshgrid(np.arange(geo_width),np.arange(geo_height))
|
|||
|
|
|
|||
|
|
lon_arr=geotrans[0]+geotrans[1]*geo_x_arr+geotrans[2]*geo_y_arr
|
|||
|
|
lat_arr=geotrans[3]+geotrans[4]*geo_x_arr+geotrans[5]*geo_y_arr
|
|||
|
|
# 行列号
|
|||
|
|
#lon_arr,lat_arr=np.meshgrid(lonlist,latlist)
|
|||
|
|
x_arr=ImageHandler.get_band_array(in_geo_x_tiff_path,1) # lon x
|
|||
|
|
y_arr=ImageHandler.get_band_array(in_geo_y_tiff_path,1) # lat y
|
|||
|
|
dem_arr=ImageHandler.get_band_array(dem_path,1)
|
|||
|
|
grid_x=np.arange(range_width)
|
|||
|
|
grid_y=np.arange(range_height)
|
|||
|
|
|
|||
|
|
x_arr=x_arr.reshape(-1).astype(np.float64)
|
|||
|
|
y_arr=y_arr.reshape(-1).astype(np.float64)
|
|||
|
|
lon_arr=lon_arr.reshape(-1).astype(np.float64)
|
|||
|
|
lat_arr=lat_arr.reshape(-1).astype(np.float64)
|
|||
|
|
dem_arr=dem_arr.reshape(-1).astype(np.float64)
|
|||
|
|
x_mask=x_arr>0
|
|||
|
|
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_table<nncols)*(rangeCoord_table>0)*(azimuthCoord_table>0)*(azimuthCoord_table<nnrows) # 注意0区域
|
|||
|
|
x=rangeCoord_table[np.where(mask==1)].reshape(-1)
|
|||
|
|
y=azimuthCoord_table[np.where(mask==1)].reshape(-1)
|
|||
|
|
z=ztd_array[np.where(mask==1)].reshape(-1)
|
|||
|
|
nx,ny=np.meshgrid(np.arange(nncols),np.arange(nnrows))
|
|||
|
|
nx=nx.reshape(-1)
|
|||
|
|
ny=ny.reshape(-1)
|
|||
|
|
ztd_range_array = scipy.interpolate.griddata((y,x),z,(ny,nx))
|
|||
|
|
ztd_range_array=ztd_range_array.reshape(nnrows,nncols)
|
|||
|
|
ImageHandler.write_imgArray(out_range_ztd_path,ztd_range_array)
|
|||
|
|
else:
|
|||
|
|
lon_table=ImageHandler.get_band_array(reflect_table_path,1)
|
|||
|
|
lat_table=ImageHandler.get_band_array(reflect_table_path,2)
|
|||
|
|
#
|
|||
|
|
im_geotrans = ImageHandler.get_dataset(ztd_array).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)
|
|||
|
|
# 斜距DLOS
|
|||
|
|
ztd_array = ImageHandler.get_band_array(in_ztd_file, 1) # 获取大气延迟表ztd
|
|||
|
|
logger.info("project ztd from geo to sar :{}".format(out_range_ztd_path))
|
|||
|
|
ztd_array=bilinearInterplor(ztd_array,sample_in_tiff_x,sample_in_tiff_y)
|
|||
|
|
ztd_array=ztd_array.reshap(shape_range)
|
|||
|
|
ImageHandler.write_imgArray(out_range_ztd_path,ztd_array)
|
|||
|
|
|
|||
|
|
|
|||
|
|
def cal_dlos(self,incidentAngle_file,in_range_ztd_path,out_Dlos_path):
|
|||
|
|
""" 计算延迟路径
|
|||
|
|
Args:
|
|||
|
|
incidentAngle_file (_type_): 入射角
|
|||
|
|
ztd_file (_type_): ztd_file
|
|||
|
|
in_lon_lat_path (_type_): in_lon_lat_path 映射表
|
|||
|
|
out_Dlos_path:
|
|||
|
|
Returns:
|
|||
|
|
_type_: _description_
|
|||
|
|
"""
|
|||
|
|
local_angle_array = ImageHandler.get_band_array(incidentAngle_file, 1) # 获取局部入射角
|
|||
|
|
ztd_array=ImageHandler.get_band_array(in_range_ztd_path, 1)
|
|||
|
|
logger.info('cal Dlos :{}'.format(out_Dlos_path))
|
|||
|
|
|
|||
|
|
dlos=ztd_array/local_angle_array
|
|||
|
|
ImageHandler.write_imgArray(out_Dlos_path,dlos)
|
|||
|
|
logger.info('DLOS cal completed !')
|
|||
|
|
|
|||
|
|
|
|||
|
|
def create_folder(folder_path):
|
|||
|
|
if not os.path.exists(folder_path):
|
|||
|
|
os.makedirs(folder_path)
|
|||
|
|
logger.info("创建文件夹：{}".format(folder_path))
|
|||
|
|
return True
|
|||
|
|
|
|||
|
|
def tropo_EAR5(inps):
|
|||
|
|
|
|||
|
|
# 坐标映射表
|
|||
|
|
reflect_lon_path=inps.lon_path #isce
|
|||
|
|
reflect_lat_path=inps.lat_path # isce
|
|||
|
|
master_pha=None
|
|||
|
|
aux_pha=None
|
|||
|
|
master_ztd_out_path = os.path.join(inps.work_spacePath, "master_ztd")
|
|||
|
|
aux_ztd_out_path = os.path.join(inps.work_spacePath, "aux_ztd")
|
|||
|
|
create_folder(master_ztd_out_path)
|
|||
|
|
create_folder(aux_ztd_out_path)
|
|||
|
|
if inps.gacos:#inps.geo2range:
|
|||
|
|
# gacos方案
|
|||
|
|
master_path=os.path.join(inps.work_spacePath,"master_ztd.tiff")
|
|||
|
|
aux_path=os.path.join(inps.work_spacePath,"aux_ztd.tiff")
|
|||
|
|
logger.info("计算主 延迟路径")
|
|||
|
|
master_ztd=Gacosdelay(master_ztd_out_path)
|
|||
|
|
master_ztd.set_waveLen(inps.wavelen)
|
|||
|
|
master_ztd.set_lon(reflect_lon_path)
|
|||
|
|
master_ztd.set_lat(reflect_lat_path)
|
|||
|
|
master_ztd.set_incidenceAngle(inps.inc_path)
|
|||
|
|
master_ztd.set_gacos_ztd(inps.mgacos,master_path)
|
|||
|
|
master_pha=master_ztd.get_delay()
|
|||
|
|
master_ztd.save()
|
|||
|
|
del master_ztd
|
|||
|
|
logger.info("计算辅 延迟路径")
|
|||
|
|
aux_ztd=Gacosdelay(aux_ztd_out_path)
|
|||
|
|
aux_ztd.set_waveLen(inps.wavelen)
|
|||
|
|
aux_ztd.set_lon(reflect_lon_path)
|
|||
|
|
aux_ztd.set_lat(reflect_lat_path)
|
|||
|
|
aux_ztd.set_incidenceAngle(inps.inc_path)
|
|||
|
|
aux_ztd.set_gacos_ztd(inps.agacos,aux_path)
|
|||
|
|
aux_pha=aux_ztd.get_delay()
|
|||
|
|
aux_ztd.save()
|
|||
|
|
del aux_ztd
|
|||
|
|
|
|||
|
|
if inps.aps:
|
|||
|
|
# 主
|
|||
|
|
logger.info("计算master 大气延迟结果")
|
|||
|
|
master_ztd=AtomZtdDelay(master_ztd_out_path)
|
|||
|
|
master_ztd.set_waveLen(inps.wavelen)
|
|||
|
|
master_ztd.get_deley_dlf(inps.MasterNCPath)
|
|||
|
|
master_ztd.get_dem(inps.DEMPath)
|
|||
|
|
master_ztd.get_lon_lat(reflect_lon_path,reflect_lat_path)
|
|||
|
|
master_ztd.set_incidenceAngle(inps.inc_path)
|
|||
|
|
master_pha=master_ztd.get_delay() # 获取
|
|||
|
|
master_ztd.save()
|
|||
|
|
del master_ztd
|
|||
|
|
|
|||
|
|
# 辅
|
|||
|
|
logger.info("计算aux 大气延迟结果")
|
|||
|
|
aux_ztd=AtomZtdDelay(aux_ztd_out_path)
|
|||
|
|
aux_ztd.set_waveLen(inps.wavelen)
|
|||
|
|
aux_ztd.get_deley_dlf(inps.AuxiliaryNCPath)
|
|||
|
|
aux_ztd.get_dem(inps.DEMPath)
|
|||
|
|
aux_ztd.get_lon_lat(reflect_lon_path,reflect_lat_path)
|
|||
|
|
aux_ztd.set_incidenceAngle(inps.inc_path)
|
|||
|
|
aux_pha=aux_ztd.get_delay()
|
|||
|
|
aux_ztd.save()
|
|||
|
|
del aux_ztd
|
|||
|
|
|
|||
|
|
# 大气延迟
|
|||
|
|
logger.info("计算大气延迟")
|
|||
|
|
out_pha_path=os.path.join(inps.outpath,'topo_pha.tiff')
|
|||
|
|
createpha(master_pha,aux_pha,inps.wavelen,out_pha_path)
|
|||
|
|
out_result_path=os.path.join(inps.outpath,os.path.basename(inps.ifg_path).replace(".tif","_ERA5.tif"))
|
|||
|
|
# 校正影像
|
|||
|
|
logger.info("校正影像")
|
|||
|
|
correct_single_ifgram(inps.ifg_path, out_pha_path,out_result_path)
|
|||
|
|
logger.info("校正结束：{}".format(out_result_path))
|
|||
|
|
pass
|