""" @Project ：microproduct @File ：AtmosphericDelayAlg.py @Function ：主程序，执行入口 @Author ：LMM @Date ：2021/09/3 14:39 @Version ：1.0.0 """ import numpy as np import netCDF4 as Nc from osgeo import gdal import math from PIL import Image import scipy.spatial import scipy.spatial.transform._rotation_groups # from scipy.interpolate import interp1d # from scipy import interpolate, integrate from scipy import interpolate, integrate # from pykrige.ok import OrdinaryKriging # import pykrige from pykrige import OrdinaryKriging from xml.etree.ElementTree import ElementTree from tool.algorithm.image.ImageHandle import ImageHandler from AtmosphericDelayAuxData import ReadImage import logging from tool.algorithm.algtools.logHandler import LogHandler import gc LogHandler.init_log_handler(r"run_log\AtmosphericDelay") logger = logging.getLogger("mylog") Rd = 287.05 # [j/kg/K] specific gas constants dry air Rv = 461.495 # [j/kg/K] water vapour k1 = 0.776*100 # [K/hPa] k2 = 0.716*100 # [K/hPa] k22 = 0.233*100 # [K/hPa] k3 = 3.75*1000*100 # [K^2/hPa] zref = 15000 # zref for integration- tropopause zincr = 15 # z increment spacing for vertical interpolation before integration vertres = 5 # vertical resolution for comparing dem to vert profiles of delay class AtmosphericDelay: """ function： 计算天顶距延迟值ZTD """ def __init__(self): pass def aps_wrf_sar(self, file, temp, re_hum, geo_height, dem_file, p_level=None): """ 计算干湿图 para: file：NC_geo_h.tif :参考数据，用来获取气象数据的经纬度、分辨率等信息 para: temp：NC_temp.tif para: re_hum：NC_re_hum.tif para: geo_height：NC_geo_h.tif para: p_level：NC_level.tif para: dem_file：dem影像数据 """ dem, max_dem = ReadImage().read_dem(dem_file) lon, lat, lon_min, lat_max, lon_max, lat_min, n_lat, n_lon, nncols, nnrows, lon_res, lat_res,\ lats, lons = ReadImage().read_lon_lat(file) lon_dem, lat_dem, lon_min_dem, lat_max_dem, lon_max_dem, lat_min_dem, n_lat_dem, n_lon_dem, nncols_dem, nnrows_dem, lon_res_dem, lat_res_dem, \ lats_dem, lons_dem = ReadImage().read_lon_lat(dem_file) min_lon_lat_dem = min(n_lon_dem, n_lat_dem) logger.info('ERA data ready !') # 1、 计算e svpw = 6.1121 * np.exp(((17.502 * (temp - 273.16)) / (240.97 + temp - 273.16))) svpi = 6.1121 * np.exp((22.587 * (temp - 273.16)) / (273.86 + temp - 273.16)) temp_bound1 = 273.16 # 0 temp_bound2 = 250.16 # -23 wgt = (temp - temp_bound2)/(temp_bound1 - temp_bound2) svp = svpi + (svpw - svpi) * (wgt ** 2) for k in range(0, 37): for i in range(0, n_lat): for j in range(0, n_lon): if temp[i, j, k].all() > 273.16: svp[i, j, k] = svpw[i, j, k] elif temp[i, j, k].all() < 250.16: svp[i, j, k] = svpi[i, j, k] else: svp[i, j, k] = svp[i, j, k] # svp = 6.11 * np.exp(2500000 / Rv * ((temp - 273.16) / (temp * 273.16))) e = re_hum*svp/100 # 2、 计算pressure pressure = p_level # {ndarray:(5,5,37)} # 3、 计算get_z g0 = 9.80665 # 重力加速度随纬度和高度而变化 # geo_height = geo_height / g0 # g = 9.80616 * (1 - 0.002637 * np.cos(2 * lats*(math.pi/180)) + 0.0000059 * (np.cos(2 * lats*(math.pi/180)))**2) # r_max = 6378137 # r_min = 6356752 # r_e = np.sqrt(1/(((np.cos(lats*(math.pi/180))**2)/(r_max**2)) + ((np.sin(lats*(math.pi/180))**2)/(r_min**2)))) # get_z = (geo_height * r_e) / (g / g0 * r_e - geo_height) get_z = geo_height/g0 cd_slices = max_dem // vertres + 1 min_lon_lat = min(n_lon, n_lat) max_lon_lat = max(n_lon, n_lat) cd_stack_dry = np.zeros((min_lon_lat, cd_slices), dtype=float) cd_stack_wet = np.zeros((min_lon_lat, cd_slices), dtype=float) # 4、 计算0-15000插值 x_i = np.array(range(0, zref, zincr)) logger.info('begin calc wet and dry of diff levels !') for m in range(0, min_lon_lat): yn = m xn = m geopo_hgt = get_z[yn, xn, :] # 位势高度 [37,] # y_p = pressure[yn, xn, :] # 压力水平level [37,] y_p = 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]) y_e = e[yn, xn, :] # 相对湿度 [37,] y_t = temp[yn, xn, :] # 温度 [37,] y_e_i = np.zeros((1, 1000), dtype=float) y_p_i = np.zeros((1, 1000), dtype=float) y_t_i = np.zeros((1, 1000), dtype=float) # 判断位势能是否在0-15000之间，避免超出范围的位势能出现插值异常 if np.min(x_i) > np.min(geopo_hgt) and np.max(x_i) < np.max(geopo_hgt): tre_yei = interpolate.interp1d(geopo_hgt, y_e) y_e_i = tre_yei(x_i) tre_yt = interpolate.interp1d(geopo_hgt, y_t) y_t_i = tre_yt(x_i) tre_yp = interpolate.interp1d(geopo_hgt, y_p) y_p_i = tre_yp(x_i) if np.min(x_i) < np.min(geopo_hgt) and np.max(x_i) < np.max(geopo_hgt): for f in range(0, 999): if x_i[f+1] > np.min(geopo_hgt) > x_i[f]: tre_yei = interpolate.interp1d(geopo_hgt, y_e) y_e_i[0, f+1:] = tre_yei(x_i[f+1:]) y_e_i[0, 0:f + 1] = y_e[np.argmin(geopo_hgt)] tre_yt = interpolate.interp1d(geopo_hgt, y_t) y_t_i[0, f+1:] = tre_yt(x_i[f+1:]) y_t_i[0, 0:f + 1] = y_t[np.argmin(geopo_hgt)] tre_yp = interpolate.interp1d(geopo_hgt, y_p) y_p_i[0, f+1:] = tre_yp(x_i[f+1:]) y_p_i[0, 0:f + 1] = y_p[np.argmin(geopo_hgt)] # {ndarray:(1:1000)} break if np.min(x_i) > np.min(geopo_hgt) and np.max(x_i) > np.max(geopo_hgt): for g in range(0, 999): if x_i[g + 1] > np.max(geopo_hgt) > x_i[g]: tre_yei = interpolate.interp1d(geopo_hgt, y_e) y_e_i[0, 0:g+1] = tre_yei(x_i[0:g+1]) y_e_i[0, g+1:] = y_e[np.argmax(geopo_hgt)] tre_yt = interpolate.interp1d(geopo_hgt, y_t) y_t_i[0, 0:g+1] = tre_yt(x_i[0:g+1]) y_t_i[0, g+1:] = y_t[np.argmax(geopo_hgt)] tre_yp = interpolate.interp1d(geopo_hgt, y_p) y_p_i[0, 0:g+1] = tre_yp(x_i[0:g+1]) y_p_i[0, g+1:] = y_p[np.argmax(geopo_hgt)] if np.min(x_i) < np.min(geopo_hgt) and np.max(x_i) > np.max(geopo_hgt): for g in range(0, 999): if x_i[g + 1] > np.min(geopo_hgt) > x_i[g]: m = g+1 y_e_i[0, 0:m] = y_e[np.argmin(geopo_hgt)] y_t_i[0, 0:m] = y_t[np.argmin(geopo_hgt)] y_p_i[0, 0:m] = y_p[np.argmin(geopo_hgt)] if x_i[g + 1] > np.max(geopo_hgt) > x_i[g]: n = g+1 y_e_i[0, n:] = y_e[np.argmax(geopo_hgt)] y_t_i[0, n:] = y_t[np.argmax(geopo_hgt)] y_p_i[0, n:] = y_p[np.argmax(geopo_hgt)] tre_yei = interpolate.interp1d(geopo_hgt, y_e) y_e_i[0, m:n] = tre_yei(x_i[m:n]) tre_yt = interpolate.interp1d(geopo_hgt, y_t) y_t_i[0, m:n] = tre_yt(x_i[m:n]) tre_yp = interpolate.interp1d(geopo_hgt, y_p) y_p_i[0, m:n] = tre_yp(x_i[m:n]) # 1、 计算湿延迟 # tmp_01 = ((k2 - (Rd * k1 / Rv)) * YeI / YTI + k3 * YeI / (YTI ** 2)) tmp_01 = k22 * y_e_i / y_t_i+k3 * y_e_i / (y_t_i ** 2) lw = np.zeros((1, 1000), dtype=float) lw1 = (10 ** (-6)) * integrate.cumtrapz(tmp_01, x_i) # (1,999) lw1 = np.fliplr(lw1) # 左右翻转(1,999) lw[0, 0:999] = lw1[0, :] # 2、 计算干延迟 tmp_02 = k1 * y_p_i / y_t_i ld = np.zeros((1, 1000), dtype=float) ld_1 = (10 ** (-6)) * integrate.cumtrapz(tmp_02, x_i) # （1,999） ld_1 = np.fliplr(ld_1) # （1,999） ld[0, 0:999] = ld_1[0, :] # ld = 2.2768 * YPI/(1-0.00266*np.cos(lats[0, m, 0]*(math.pi/180))-0.00028*XI) # ld计算方法2 横着（1000，） # ld =(10 ** (-6)) *(k1*Rd/g0)*(YPI-YPI[0,999]) # ld计算方法3 # 3、在(0,maxdem)上做插值 cdi = np.array(range(0, max_dem, vertres)) xi_range = x_i ld_range = ld[0, :] lw_range = lw[0, :] # ld_i = interpolate.interp1d(xi_range, ld_range, kind="cubic") ld_i = interpolate.interp1d(xi_range, ld_range) # lw_i = interpolate.interp1d(xi_range, lw_range, kind="cubic") lw_i = interpolate.interp1d(xi_range, lw_range) ldi = np.zeros((max_dem//vertres+1, 1), dtype=float) lwi = np.zeros((max_dem//vertres+1, 1), dtype=float) for i in range(0, max_dem // vertres): ldi[i, 0] = ld_i(cdi[i]) lwi[i, 0] = lw_i(cdi[i]) cd_stack_dry[m, :] = ldi[:, 0].T cd_stack_wet[m, :] = lwi[:, 0].T logger.info('finish calc wet and dry delay of diff levels !') cdstack_interp_dry_wet = np.zeros((min_lon_lat_dem, min_lon_lat_dem, cd_stack_dry.shape[1]), dtype=float) # 修改插值网格分辨率 # cdstack_interp_dry_wet = np.zeros((min_lon_lat, min_lon_lat, cd_stack_dry.shape[1]), dtype=float) # np.zeros((max_lon_lat, 1), dtype=float) # 创建data01、data02 data01 = np.zeros((min_lon_lat, 1), dtype=float) data02 = np.zeros((min_lon_lat, 1), dtype=float) new_lon = np.zeros((n_lon, 1), dtype=float) new_lat = np.zeros((n_lat, 1), dtype=float) new_lon[:, 0] = lon new_lat[:, 0] = lat data01[0:min_lon_lat, 0] = new_lon[0:min_lon_lat, 0] data02[0:min_lon_lat, 0] = new_lat[0:min_lon_lat, 0] data01 = data01.reshape(data01.shape[0], 1) # 水平转垂直 data02 = data02.reshape(data02.shape[0], 1) # 水平转垂直 lat_max = lat_min+(min_lon_lat-1) * (-lat_res) lon_max = lon_min+(min_lon_lat-1) * lon_res logger.info('start interpolation dry and wet delay...') # 在网格上进行克里金插值 # 1. 先插值到 气象数据给的分辨率 try: for i in range(0, cd_stack_dry.shape[1] - 1): data03 = cd_stack_dry[:, i] + cd_stack_wet[:, i] cdstack_interp_dry_wet[:, :, i] = self.kringing(data01, data02, data03, lon_min, lat_max, lon_max, lat_min, lon_res_dem, -lat_res_dem, min_lon_lat_dem) except: logger.info('calc kringing error !') logger.info('start calc rbf ... !') for i in range(0, cd_stack_dry.shape[1] - 1): data03 = cd_stack_dry[:, i] + cd_stack_wet[:, i] cdstack_interp_dry_wet[:, :, i] = self.calc_rbf(data01, data02, data03, lon_min, lat_max, lon_max, lat_min, lon_res_dem, -lat_res_dem, min_lon_lat_dem) logger.info('calc kringing interv finish !') return cdstack_interp_dry_wet # try: # for i in range(0, cd_stack_dry.shape[1]-1): # data03 = cd_stack_dry[:, i]+cd_stack_wet[:, i] # cdstack_interp_dry_wet[:, :, i] = self.kringing(data01, data02, data03, lon_min, lat_max, lon_max, # lat_min, lon_res, -lat_res, min_lon_lat) # except: # logger.info('calc kringing error !') # logger.info('start calc rbf ... !') # for i in range(0, cd_stack_dry.shape[1]-1): # data03 = cd_stack_dry[:, i]+cd_stack_wet[:, i] # cdstack_interp_dry_wet[:, :, i] = self.calc_rbf(data01, data02, data03, lon_min, lat_max, lon_max, # lat_min, lon_res, -lat_res, min_lon_lat) # logger.info('calc kringing interv finish !') @staticmethod def calc_ztd(dry_wet_file_name, dem_file): """ 计算ztd wet_file_name:干湿校正文件路径 dem_file: dem文件路径 """ cdstack_interp_dry_wet = ImageHandler().get_all_band_array(dry_wet_file_name) dem = ImageHandler().get_all_band_array(dem_file) nncols = ImageHandler().get_img_width(dem_file) nnrows = ImageHandler().get_img_height(dem_file) wet_dry_wid = ImageHandler().get_img_width(dry_wet_file_name) dem_wid = ImageHandler().get_img_width(dem_file) wet_dry_hgt = ImageHandler().get_img_height(dry_wet_file_name) dem_hgt = ImageHandler().get_img_height(dem_file) wetcorrection = np.zeros((nnrows, nncols), dtype=float) rounddem = np.round(dem/vertres) rounddem[np.where(dem < 0)] = 0 ratio_wid = (dem_wid/wet_dry_wid) ratio_hgt = (dem_hgt/wet_dry_hgt) for i in range(1, wet_dry_hgt+1): for j in range(1, wet_dry_wid+1): m1 = round((i - 1) * ratio_hgt) m2 = round(i * ratio_hgt) n1 = round((j - 1) * ratio_wid) n2 = round(j * ratio_wid) t1_array02 = np.zeros((m2-m1, n2-n1), dtype=float) for h in range(0, cdstack_interp_dry_wet.shape[2]): t1_array = rounddem[m1:m2, n1:n2] t1_array01 = np.where(t1_array == h, 1, 0) t1_array01 = np.where(t1_array01 > 0, cdstack_interp_dry_wet[i-1, j-1, h], 0) t1_array02 = t1_array02+t1_array01 wet_block_array = t1_array02 wetcorrection[m1:m2, n1:n2] = wet_block_array ztd = wetcorrection logger.info('calc ztd finish !') return ztd @staticmethod def read_nc(infile): """ 读取气象数据 infile：气象数据路径 """ nc_data_obj = Nc.Dataset(infile) keys = nc_data_obj.variables.keys() if 'longitude' in keys: lons = nc_data_obj.variables['longitude'][:] lats = nc_data_obj.variables['latitude'][:] elif 'lon' in keys: lons = nc_data_obj.variables['lon'][:] lats = nc_data_obj.variables['lat'][:] else: logger.info("Specify variable names in the .nc file !") temp = np.array(nc_data_obj.variables["t"]) # temp in K 温度 # temp =(temp-255.6357)/0.00153279 hum = np.array(nc_data_obj.variables["r"]) # relative humidity in percent 湿度 # hum = (hum-41.8782447)/0.0012780 geo_h = np.array(nc_data_obj.variables["z"]) # Geopotential Height in m 重力势 # geo_h =(H-237742.3474)/7.243160 p_levs = np.array(nc_data_obj.variables["level"]) nc_data_obj.close() # 关闭文件 return lons, lats, temp, hum, geo_h, p_levs @staticmethod def read_dem(demfile): """ 读取dem数据 demfile：dem数据路径 """ dem = np.array(Image.open(demfile)) maxdem = math.ceil(np.max(dem)/100)*100+50 # 四舍五入向上取值 return dem, maxdem @staticmethod def copy_channel(num_lat, num_lon, air_temp, rel_hum, h, lons, lats, plevs): """ 改变数据维度 """ temp_1 = np.zeros((num_lat, num_lon, 37), dtype=float) hum_1 = np.zeros((num_lat, num_lon, 37), dtype=float) h_1 = np.zeros((num_lat, num_lon, 37), dtype=float) lat_1 = np.zeros((1, num_lat, 37), dtype=float) lon_1 = np.zeros((1, num_lon, 37), dtype=float) p_level_1 = np.zeros((num_lat, num_lon, 37), dtype=float) for i in range(0, 37): a = air_temp[0, i, :, :] temp_1[0:num_lat, 0:num_lon, i] = a b = rel_hum[0, i, :, :] hum_1[0:num_lat, 0:num_lon, i] = b c = h[0, i, :, :] h_1[:, :, i] = c del air_temp, rel_hum, h air_temp, rel_hum, h = temp_1, hum_1, h_1 for i in range(0, 37): d = lons lon_1[0, :, i] = d e = lats lat_1[0, :, i] = e del lons, lats lons, lats = lon_1, lat_1 # Plevs = Plevs/100 for i in range(0, 37): p_level_1[:, :, i] = p_level_1[:, :, i] + plevs[i] # print("输出Plevs1", Plevs1[:, :, i]) p_level = p_level_1 return air_temp, rel_hum, h, lons, lats, p_level @staticmethod def calc_rbf(data01, data002, data003, lon_min, lat_max, lon_max, lat_min, lon_res, lat_res, min_lon_lat): """ rbf插值 """ # gridx = np.arange(lon_min, lon_max, lon_res) # 三个参数的意思：范围0.0 - 0.6 ，每隔0.1划分一个网格 # gridy = np.arange(lat_min, lat_max, lat_res) gridx = np.linspace(lon_min, lon_max, min_lon_lat) # 三个参数的意思：范围(lon_min, lon_max) ,每隔Lon_Res划分一个网格 gridy = np.linspace(lat_min, lat_max, min_lon_lat) func = interpolate.Rbf(data01, data002, data003, function='multiquadric') znew = func(gridx, gridy) return znew @staticmethod def kringing(data01, data02, data03, lon_min, lat_max, lon_max, lat_min, lon_res, lat_res, min_lon_lat): """ 克里金插值方法 """ # gridx = np.arange(lon_min, lon_max+lon_res, lon_res) # 三个参数的意思：范围(lon_min, lon_max) ,每隔Lon_Res划分一个网格 # gridy = np.arange(lat_min, lat_max+lat_res, lat_res) # ok3d = ok.OrdinaryKriging(data01, data02, data03, variogram_model="linear") # 模型 # ok3d = ok.OrdinaryKriging(data01, data02, data03) # 模型 gridx = np.linspace(lon_min, lon_max, min_lon_lat) # 三个参数的意思：范围(lon_min, lon_max) ,每隔Lon_Res划分一个网格 gridy = np.linspace(lat_min, lat_max, min_lon_lat) # ok3d = pykrige.OrdinaryKriging(data01, data02, data03) # 模型 ok3d = OrdinaryKriging(data01, data02, data03) k3d1, ss3d = ok3d.execute("grid", gridx, gridy) # k3d1是结果，给出了每个网格点处对应的值 # print(np.round(k3d1, 6)) # 小数点后6位 return k3d1 @staticmethod def write_dry_wet_tif(ref_image_dem, ref_image_nc, in_array, out_path): """ 存为tif数据 ref_image：参考影像路径 in_array：输入数组 outpath：输出路径 """ dataset_dem = gdal.Open(ref_image_dem) # 获取dem宽高及分辨率 x_size = dataset_dem.RasterXSize y_size = dataset_dem.RasterYSize trans_dem = dataset_dem.GetGeoTransform() dataset_nc = gdal.Open(ref_image_nc) # 获取气象数据坐标信息 trans_nc = dataset_nc.GetGeoTransform() tran_geo = [trans_nc[0], trans_dem[1], 0, trans_nc[3], 0, trans_dem[5]] if len(in_array.shape) == 3: im_height, im_width, im_bands = in_array.shape else: im_bands, (im_height, im_width) = 1, in_array.shape pro = dataset_nc.GetProjection() gdal.AllRegister() drive = gdal.GetDriverByName("GTiff") newfile = drive.Create(out_path, y_size, y_size, im_bands, gdal.GDT_Float32) newfile.SetProjection(pro) newfile.SetGeoTransform(tran_geo) if im_bands == 1: newfile.GetRasterBand(1).WriteArray(in_array) # 写入数组数据 else: for i in range(0, im_bands): newfile.GetRasterBand(i + 1).WriteArray(in_array[:, :, i]) # dataset.GetRasterBand(i + 1).WriteArray(im_data[i]) del newfile @staticmethod def write_tif(ref_image, in_array, out_path): """ 存为tif数据 ref_image：参考影像路径 in_array：输入数组 outpath：输出路径 """ dataset = gdal.Open(ref_image) x_size = dataset.RasterXSize y_size = dataset.RasterYSize if len(in_array.shape) == 3: im_height, im_width, im_bands = in_array.shape else: im_bands, (im_height, im_width) = 1, in_array.shape pro = dataset.GetProjection() geo = dataset.GetGeoTransform() gdal.AllRegister() drive = gdal.GetDriverByName("GTiff") newfile = drive.Create(out_path, x_size, y_size, im_bands, gdal.GDT_Float32) newfile.SetProjection(pro) newfile.SetGeoTransform(geo) if im_bands == 1: newfile.GetRasterBand(1).WriteArray(in_array) # 写入数组数据 else: for i in range(0, im_bands): newfile.GetRasterBand(i + 1).WriteArray(in_array[:, :, i]) # dataset.GetRasterBand(i + 1).WriteArray(im_data[i]) del newfile @staticmethod def write_ztd_tif(ref_image_dem, ref_image_nc, in_array, out_path): """ 存为tif数据 ref_image：参考影像路径 in_array：输入数组 outpath：输出路径 """ dataset_dem = gdal.Open(ref_image_dem) x_size = dataset_dem.RasterXSize y_size = dataset_dem.RasterYSize trans_dem = dataset_dem.GetGeoTransform() dataset_nc = gdal.Open(ref_image_nc) # 获取气象数据坐标信息 trans_nc = dataset_nc.GetGeoTransform() tran_geo = [trans_nc[0], trans_dem[1], 0, trans_nc[3], 0, trans_dem[5]] if len(in_array.shape) == 3: im_height, im_width, im_bands = in_array.shape else: im_bands, (im_height, im_width) = 1, in_array.shape pro = dataset_nc.GetProjection() gdal.AllRegister() drive = gdal.GetDriverByName("GTiff") newfile = drive.Create(out_path, x_size, y_size, im_bands, gdal.GDT_Float32) newfile.SetProjection(pro) newfile.SetGeoTransform(tran_geo) if im_bands == 1: newfile.GetRasterBand(1).WriteArray(in_array) # 写入数组数据 else: for i in range(0, im_bands): newfile.GetRasterBand(i + 1).WriteArray(in_array[:, :, i]) # dataset.GetRasterBand(i + 1).WriteArray(im_data[i]) del newfile @staticmethod def cal_phase_bias(local_file, ztd_file, waveLength): """ function: 计算延迟转换,输出像素-DLOS表 param: local_file : 输入-局部入射角文件路径 param: ztd_file : 输入-大气延迟表 ZTD路径 param: waveLength : 影像对应的波长 """ ztd_array = ImageHandler.get_band_array(ztd_file, 1) # 获取大气延迟表ztd local_angle_array = ImageHandler.get_band_array(local_file, 1) # 获取局部入射角 dlos = (1 / np.cos(local_angle_array)) * ztd_array # 计算公式 phase_bias = 4 * math.pi /waveLength * dlos # 计算主影像相位误差 return dlos#phase_bias def complex_phase_removal(self, PhaseError_m, PhaseError_a, intervene_sar,waveLength): """ function: isce输出的复数干涉图 除以e^(0+wj),得到去除延迟相位的复数干涉图 param: PhaseError_m: 主影像相位误差图 param: PhaseError_a：辅影像相位误差图 param: intervene_sar：原始干涉图 """ phase_bias_m = ImageHandler.get_band_array(PhaseError_m, 1) # 计算主影像相位误差 phase_bias_a = ImageHandler.get_band_array(PhaseError_a, 1) # 计算辅影像相位误差 bias =(4 * math.pi /waveLength)*(phase_bias_m-phase_bias_a) # 延迟相位 sar_phase = ImageHandler.get_band_array(intervene_sar, 1) complex64_sar_phase = sar_phase.astype(np.complex64) complex_phase = bias * (1j) fl64_com = complex_phase.astype(np.complex64) # 另存为complex64 bias_result = complex64_sar_phase/np.exp(fl64_com) bias_result[np.isnan(bias_result)] = 0 + 0j return bias_result @staticmethod def get_radarwavelength(IW1): """ 读取波长 :param IW1: :return: """ tree = ElementTree() tree.parse(IW1) root = tree.getroot() RadarCenterFrequency = root.find('component').find('component').find('component') element_trees = list(RadarCenterFrequency) value = 0 for element in element_trees: if len(element.attrib) == 1: if element.attrib["name"]=="radarwavelength": value=float(element.find('value').text) return value @staticmethod def get_polarization(IW1): """ 读取极化方式 :param IW1: :return: """ tree = ElementTree() tree.parse(IW1) root = tree.getroot() RadarCenterFrequency = root.find('component').find('component').find('component') element_trees = list(RadarCenterFrequency) value = 0 for element in element_trees: if len(element.attrib) == 1: if element.attrib["name"]=="polarization": value=element.find('value').text return value