GACOS大气改正python实现
该代码共有8个部分。分别是数据读取、头文件读取、ztd数据裁剪、趋势项去除、相位包裹、make correction、ztd转los以及主程序
我为什么要写这个代码呢？原因有二：
首先，因为我个人不喜欢使用matlab，用着不顺手。
其次，我想锻炼一下自己的编程能力，所以就着手了并完成了这样一个问题。
注：代码有任何问题，请私信我，接受建议，谢谢！！！
注：转载请标明出处，谢谢！！！
1.数据读取
常规的二进制文件读取。可以自己先调试数据是否读取正确，可视化一下，看是否正常。代码命名为read_binary.py
import numpy as np
import struct
def xshow(filename, nx, nz):
f = open(filename, 'rb')
pic = np.zeros((nx, nz))
for i in range(nx):
for j in range(nz):
data = f.read(4)
elem = struct.unpack("f", data)[0]
pic[i][j] = elem
f.close()
return pic
2.ztd头文件读取
读取的就是下载的gacos头文件数据（.rsc）。代码命名为read_rsc.py
def read_rsc(file):
file_path = file
with open(file_path) as f:
count = len(open(file_path, 'r').readlines())
for i in range(count):
line = f.readline().strip()
line = line.split()
if line[0] == 'WIDTH':
width = line[1]
width = int(width)
elif line[0] == 'FILE_LENGTH':
file_length = line[1]
file_length = int(file_length)
elif line[0] == 'XMIN':
xmin = line[1]
xmin = int(xmin)
elif line[0] == 'XMAX':
xmax = line[1]
xmax = int(xmax)
elif line[0] == 'YMIN':
ymin = line[1]
ymin = int(ymin)
elif line[0] == 'YMAX':
ymax = line[1]
ymax = int(ymax)
elif line[0] == 'X_FIRST':
x_first = line[1]
x_first = float(x_first)
elif line[0] == 'Y_FIRST':
y_first = line[1]
y_first = float(y_first)
elif line[0] == 'X_STEP':
x_step = line[1]
x_step = float(x_step)
elif line[0] == 'Y_STEP':
y_step = line[1]
y_step = float(y_step)
else:
pass
return width, file_length, xmin, xmax, ymin, ymax, x_first, y_first, x_step, y_step
3.ztd数据裁剪
由于ztd数据的分布范围更广，需要将其裁剪到与InSAR数据相同的范围，因此需要进行裁剪。代码命名为cut_image.py
import struct
import matplotlib.pyplot as plt
from read_binary import xshow
import numpy as np
import pandas as pd
# 根据rsc头文件进行裁剪
def read_rsc(file):
file_path = file
with open(file_path) as f:
count = len(open(file_path, 'r').readlines())
# print(count)
for i in range(count):
line = f.readline().strip()
line = line.split()
if line[0] == 'WIDTH':
width = line[1]
width = int(width)
# print('width is {}'.format(width))
elif line[0] == 'FILE_LENGTH':
length = line[1]
length = int(length)
# print('file_length is {}'.format(length))
elif line[0] == 'X_FIRST':
x_first = line[1]
x_first = float(x_first)
# print('x_first is {}'.format(x_first))
elif line[0] == 'Y_FIRST':
y_first = line[1]
y_first = float(y_first)
# print('y_first is {}'.format(y_first))
elif line[0] == 'X_STEP':
x_step = line[1]
x_step = float(x_step)
# print('x_step is {}'.format(x_step))
elif line[0] == 'Y_STEP':
y_step = line[1]
y_step = float(y_step)
# print('y_step is {}'.format(y_step))
else:
pass
return width, length, x_first, y_first, x_step, y_step
def cut_image(dir, out_path, filename, maxlat, len_l, minlon, wid_l, outfilename):
file_path = dir + '\\' + filename
rsc_file_path = dir + '\\' + filename + '.rsc'
cut_outpath = out_path + '\\' + filename + outfilename
width, length, x_first, y_first, x_step, y_step = read_rsc(rsc_file_path)
# read ztd data
ztd_data = xshow(file_path, length, width)
# 可视化
plt.imshow(ztd_data)
plt.show()
# generate none lat and lon(same as ztd)
lat = np.zeros((length, width))
lon = np.zeros((length, width))
# 使用循环将空lat、lon矩阵按照间隔写入数据
for i in range(length):
for j in range(width):
lat[i, j] = y_first + y_step * i
lon[i, j] = x_first + x_step * j
index_maxlat = round((maxlat - y_first) / y_step)
index_minlat = index_maxlat + len_l - 1
index_minlon = round((minlon - x_first) / x_step)
index_maxlon = index_minlon + wid_l - 1
# 检查裁剪的数据是否正确
if index_maxlat < 1 or index_minlat > length:
raise ValueError('cut_image: the input image is smaller than the output!')
if index_minlon < 1 or index_maxlon > width:
raise ValueError('cut_image: the input image is smaller than the output!')
ztd_data_cut = ztd_data[index_maxlat:index_minlat + 1, index_minlon: index_maxlon + 1]
# 可视化裁剪后数据
plt.imshow(ztd_data_cut)
plt.show()
# 将裁剪后数据写入文件
cut_data = pd.DataFrame(data=None, columns=['cut'])
cut_data['cut'] = ztd_data_cut.flatten()
# cut_data.to_csv(cut_outpath + '.csv', index=False)
with open(cut_outpath, 'wb') as f:
for i in cut_data['cut'].values.tolist():
a = struct.pack('f', i)
f.write(a)
# f.close()
print("cut file Done")
return ztd_data_cut
4.趋势项去除
这里我理解的就是轨道相位。虽然我这里是按照gacos官方复现的，但是不是特别理解他的去除方法原理。
我对于去除轨道趋势项的理解是在一般情况下利用距离向、方位向与相位建模，求出他们的函数关系，再去除。代码命名为remove_planar.py
"""
Remove planar trend of an interferogram
The planar is defined as:
Trend = a0 + a1 * X + a2 * Y
Parameters are regressed using all available phases by the above equation
and removed from the phases.
"""
import numpy as np
import matplotlib.pyplot as plt
def remove_planar(corrected):
num = 0
sumx = 0
sumy = 0
sumxy = 0
sumx2 = 0
sumy2 = 0
sumsx = 0
sumsy = 0
sums = 0
# data = np.array(corrected)
data = corrected.flatten() # 为了获取index，所以将数据展平，后续操作还是使用corrected
# data = corrected
zero_index = np.where(data == 0)
len, wid = np.shape(corrected)
for i in range(len):
for j in range(wid):
if corrected[i][j] == 0 or np.isnan(corrected[i][j]): # 二维数组访问元素data[i][j]第i行第j列,使用continue语句剔除nan与0
continue
gridx = j
gridy = i
phase = corrected[i][j]
sumx = sumx + gridx
sumy = sumy + gridy
sumx2 = sumx2 + gridx * gridx
sumy2 = sumy2 + gridy * gridy
sumxy = sumxy + gridx * gridy
sumsx = sumsx + phase * gridx
sumsy = sumsy + phase * gridy
sums = sums + phase
num = num + 1
A = np.zeros((3, 3))
A[0, 0] = sumx2
A[0, 1] = sumxy
A[0, 2] = sumx
A[1, 0] = sumxy
A[1, 1] = sumy2
A[1, 2] = sumy
A[2, 0] = sumx
A[2, 1] = sumy
A[2, 2] = num
L = np.zeros(3)
L[0] = sumsx
L[1] = sumsy
L[2] = sums
L = np.transpose(L)
R = np.dot(np.linalg.inv(A), L)
tdata = np.zeros((len, wid))
for i in range(len):
for j in range(wid):
if corrected[i][j] == 0:
pass
gridx = j
gridy = i
tdata[i][j] = corrected[i][j] - R[0]*gridx - R[1] * gridy -R[2]
new_tdata = tdata.flatten()
new_tdata[zero_index] = 0
return new_tdata
5.相位包裹
同上，我也是按照gacos官方复现的。我唯一的感觉就是这相位包裹是否过于简单？？代码命名为wrap_matrix.py
"""
Wrap the input matrix into [vmin,vmax]
vmin,vmax不同，wrap的程度不同,设置为[-3.14,3.14]
"""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
def wrap_matrix(data, vmin, vmax):
# 获取nan值的索引
index_data = np.array(np.where(np.isnan(data.flatten())))
# index_data = phase_data[phase_data['phase'].isna()].index
# 获取小于vmin值的索引
min_index = np.argwhere(data.flatten() <= vmin)
# min_index = phase_data[phase_data['phase'] < vmin].index
size = min_index.size
output = data
output_flatten = output.flatten()
while min_index.size != 0:
# 数据均是数组--根据min_index的位置将其数值处理--不包括nan
output_flatten[min_index] = output_flatten[min_index] + (vmax-vmin)
min_index = np.argwhere(output_flatten <= vmin)
max_index = np.argwhere(data.flatten() > vmax)
while max_index.size != 0:
output_flatten[max_index] = output_flatten[max_index] - (vmax-vmin)
max_index = np.argwhere(output_flatten > vmax)
# output.flatten()[index_data] = np.nan
plt.imshow(output_flatten.reshape(1200, 1200))
plt.show()
return output_flatten
6.make correction
使用gacos进行改正。代码命名为make_correction.py
"""
实现参考点功能
思路：
1.获得数组中nan,0值的索引
2.将0值转为nan
3.计算所有非nan的均值，
4.将原始相位的绝对值减去均值，找到其最小的那个值，并得到其索引
4.根据索引获得其相位值
5.将相位每一个值减去该最小值
6.dztd同上一步
"""
import struct
import numpy as np
import pandas as pd
from read_binary import xshow
import matplotlib.pyplot as plt
import os
import cut_image
import read_rsc
import remove_planar
import wrap_matrix
import math
import operator
def make_correction(dir, procecss_dir, save_dir, ifg_name, ztd1_name, ztd2_name, ref_lat, ref_lon, wavelength,
unit_phase, unit_out, fln_elev, elev, isplanar, iswrap):
# -------------------------------------------------get some parameters----------------------------------------------
ifg_header = dir + '\\' + ifg_name + '.rsc'
width, length, xmin, xmax, ymin, ymax, x_first, y_first, x_step, y_step = read_rsc.read_rsc(ifg_header)
# -------------------cut image to be exactly the same with interferogram，data format:ndarray-----------------------
cut_ztd1 = cut_image.cut_image(dir=dir, out_path=procecss_dir, filename=ztd1_name, maxlat=y_first, len_l=length, minlon=x_first, wid_l=width, outfilename='1')
cut_ztd2 = cut_image.cut_image(dir=dir, out_path=procecss_dir, filename=ztd2_name, maxlat=y_first, len_l=length, minlon=x_first, wid_l=width, outfilename='2')
if os.path.exists(fln_elev):
cut_elev = cut_image.cut_image(dir=procecss_dir, out_path=procecss_dir, filename=fln_elev, maxlat=y_first, len_l=length, minlon=x_first, wid_l=width, outfilename='elev')
# ---------------------------------------------read phase and convert to meters-------------------------------------
phase_data = xshow(filename=dir + '\\' + ifg_name, nx=width, nz=length) # ndarray[cols, lines]
# max_d = max(phase_data[~np.isnan(phase_data)])
if unit_phase == 'p':
ifg_data = phase_data * wavelength / (4 * np.pi)
max_da = max(ifg_data[~np.isnan(ifg_data)])
elif unit_phase == 'c':
ifg_data = phase_data / 100
else:
raise ValueError('please check the read unit')
# -------------------------------------------read ztd files from GACOS----------------------------------------------
ztd1_data = cut_ztd1 # data type：ndarray
ztd2_data = cut_ztd2 # data type：ndarray
# ------------------------------------------------difference ztd----------------------------------------------------
diff_ztd = ztd2_data - ztd1_data # data type：ndarray
# 将数据转为dataframe，便于后续处理
ifg_dataframe = pd.DataFrame(data=None, columns=['ifg'])
diff_ztd_dataframe = pd.DataFrame(data=None, columns=['ztd'])
ifg_dataframe['ifg'] = ifg_data.flatten()
diff_ztd_dataframe['ztd'] = diff_ztd.flatten()
# ------------------------------先将0值转为nan值（应该是为了参考点的设置，需要去除0），然后再去除nan值--------------------------
ifg_dataframe['ifg'].replace(0, np.nan, inplace=True)
ifg_new_dataframe = ifg_dataframe.dropna()
# ------------------------------------------------set reference point-----------------------------------------------
min_index = 0
ifg_array = np.array(ifg_new_dataframe['ifg'])
diff_ztd_array = np.array(diff_ztd_dataframe['ztd'])
if ref_lon == 0 and ref_lat == 0:
mean = np.mean(ifg_array)
# print(mean)
# sub_array = np.min(np.abs(ifg_array) - mean) # 这里的绝对值很重要！
sub_array = np.abs(ifg_array - mean)
# sub_array = np.abs(ifg_dataframe['ifg'].values - mean)
# 将sub_array写入原始数据中，方便找到对应的index
sub_ifg_dataframe = pd.DataFrame(pd.Series(index=ifg_dataframe[~ifg_dataframe['ifg'].isna()].index, data=sub_array))
sub_ifg_dataframe['sub_ifg'] = sub_array.flatten()
sub_new_array = sub_ifg_dataframe['sub_ifg'].values
min_ztd_index = sub_ifg_dataframe['sub_ifg'].idxmin()
min_phase_index, min_number = min(enumerate(sub_new_array), key=operator.itemgetter(1))
# min_index = np.argmin(sub_array)
# max__d = max(sub_array)
# print(min_index)
else:
pass
change_ifg_array = ifg_array - ifg_array[min_phase_index] # 整幅干涉图中所有相位值减去参考点相位，ifg_array维度是一维，可以这样取值
min_phase = ifg_array[min_phase_index]
# max_ddd = max(change_ifg_array)
# mean_dd = np.mean(change_ifg_array)
# ddd = ifg_array[min_index]
# print(change_ifg_array)
dztd_data = diff_ztd_array - diff_ztd_array[min_ztd_index]
min_ztd = diff_ztd_array[min_ztd_index]
# ------------------------------------------------设置输出的单位-------------------------------------------------------
if unit_out == 'p':
new_ifg_array = change_ifg_array / (wavelength / (4 * np.pi))
# max_dddd = max(new_ifg_array)
new_dztd_data = dztd_data / (wavelength / (4 * np.pi))
elif unit_out == 'c':
new_ifg_array = change_ifg_array * 100.0
new_dztd_data = dztd_data * 100.0
else:
raise ValueError('please check the output unit')
# ----------------------------------------------read elevation data-------------------------------------------------
if os.path.exists(fln_elev):
elevdata = xshow(filename=dir + fln_elev + 'elev', nx=width, nz=length)
dztd = np.divide(new_dztd_data, np.sin(elevdata))
else:
elev = elev * 3.14159265 / 180.0
dztd = np.divide(new_dztd_data, np.sin(elev))
dztd = dztd.reshape(width, length)
# -----------------------------------将转换单位的dztd转为dataframe,便于后续保存------------------------------------------
dztd_dataframe = pd.DataFrame(data=None, columns=['dztd'])
dztd_dataframe['dztd'] = dztd.flatten()
# -----------------------------------将转换单位的ifg转成dataframe--便于与dztd作差----------------------------------------
new_ifg_dataframe = pd.DataFrame(pd.Series(index=ifg_dataframe[~ifg_dataframe['ifg'].isna()].index, data=new_ifg_array))
new_ifg_dataframe['new_ifg'] = new_ifg_array.flatten()
# max_ddd = max(new_ifg_dataframe['new_ifg'])
print()
# ------------------------------------------------correction--------------------------------------------------------
ifg_dataframe['corrected'] = new_ifg_dataframe['new_ifg'] - dztd_dataframe['dztd']
ifg_dataframe['phase'] = new_ifg_dataframe['new_ifg']
# --------------------------------------------------保存ztd解缠相位------------------------------------------------------
with open(save_dir + '\\' + ifg_name + '-phase-dztd', 'wb') as f:
for i in ifg_dataframe['corrected'].values.tolist():
a = struct.pack('f', i)
f.write(a)
print('-------------------')
print("corrected file Done")
print('-------------------')
# 根据原始相位将dztd映射到相同的位置
phase_label = ifg_dataframe['corrected'].values
dztd_label = dztd_dataframe['dztd'].values
for i in range(len(phase_label)):
if math.isnan(phase_label[i]):
dztd_label[i] = np.nan
with open(save_dir + '\\' + ifg_name +'-dztd', 'wb') as f:
for i in dztd_dataframe['dztd'].values.tolist():
a = struct.pack('f', i)
f.write(a)
print('-------------------')
print("dztd file Done")
print('-------------------')
with open(save_dir + '\\' + ifg_name, 'wb') as f:
for i in ifg_dataframe['phase'].values.tolist():
a = struct.pack('f', i)
f.write(a)
print('-------------------')
print("phase file Done")
print('-------------------')
# -------------------------------------------------是否去除趋势项------------------------------------------------------
if isplanar != 0:
print('-------------------')
print('开始去除趋势项')
print('-------------------')
corrected_array = np.array(ifg_dataframe['corrected']).reshape(width, length)
detrend = remove_planar.remove_planar(corrected_array)
ifg_dataframe['detrend'] = detrend.flatten()
with open(save_dir + '\\' + ifg_name + '-phase-dztd-planar', 'wb') as f:
for i in ifg_dataframe['detrend'].values.tolist():
a = struct.pack('f', i)
f.write(a)
print('-------------------')
print("detrend file Done")
print('-------------------')
# --------------------------------------------------可视化-----------------------------------------------------------
print('-------------------')
print('可视化GACOS改正结果')
print('-------------------')
fig, ax = plt.subplots(2, 2)
ax1 = ax[0][0].imshow(np.array(ifg_dataframe['phase']).reshape(width, length), cmap='jet')
cb1 = plt.colorbar(ax1, ax=ax[0][0])
ax2 = ax[0][1].imshow(np.array(dztd_dataframe['dztd']).reshape(width, length), cmap='jet')
cb2 = plt.colorbar(ax2, ax=ax[0][1])
ax3 = ax[1][0].imshow(np.array(ifg_dataframe['corrected']).reshape(width, length), cmap='jet')
cb3 = plt.colorbar(ax3, ax=ax[1][0])
ax4 = ax[1][1].imshow(np.array(ifg_dataframe['detrend']).reshape(width, length), cmap='jet')
cb4 = plt.colorbar(ax4, ax=ax[1][1])
cb1.ax.set_title('(rad)')
cb2.ax.set_title('(rad)')
cb3.ax.set_title('(rad)')
cb4.ax.set_title('(rad)')
fig.suptitle('GACOS corrected results')
ax[0][0].set_title('phase')
ax[0][1].set_title('dztd')
ax[1][0].set_title('phase-dztd')
ax[1][1].set_title('phase-dztd-planar')
plt.show()
# -------------------------------------------------------保存包裹相位-------------------------------------------------
if iswrap > 0:
print('-------------------')
print('开始将相位包裹')
print('-------------------')
dztd_wrap = wrap_matrix.wrap_matrix(np.array(dztd_dataframe['dztd'].values).reshape(width, length),
-iswrap, iswrap)
phase_wrap = wrap_matrix.wrap_matrix(np.array(ifg_dataframe['ifg'].values).reshape(width, length), -iswrap,
iswrap)
if isplanar != 0:
corr_wrap = wrap_matrix.wrap_matrix(np.array(ifg_dataframe['detrend'].values).reshape(width, length),
-iswrap, iswrap)
with open(save_dir + '\\' + ifg_name + '-phase-dztd-planar_wrap', 'wb') as f:
for i in corr_wrap:
a = struct.pack('f', i)
f.write(a)
print('-------------------')
print("detrend_wrap file Done")
print('-------------------')
else:
corr_wrap = wrap_matrix.wrap_matrix(np.array(ifg_dataframe['corrected'].values).reshape(width, length),
-iswrap, iswrap)
with open(save_dir + '\\' + ifg_name + '-phase-dztd_wrap', 'wb') as f:
for i in corr_wrap:
a = struct.pack('f', i)
f.write(a)
print('-------------------')
print("phase-dztd_wrap file Done")
print('-------------------')
with open(save_dir + '\\' + ifg_name + '-dztd_wrap', 'wb') as f:
for i in dztd_wrap:
a = struct.pack('f', i)
f.write(a)
print('-------------------')
print("dztd_wrap file Done")
print('-------------------')
with open(save_dir + '\\' + ifg_name + '-phase_wrap', 'wb') as f:
for i in phase_wrap:
a = struct.pack('f', i)
f.write(a)
print('-------------------')
print("phased_wrap file Done")
print('-------------------')
7.主程序
实现单干涉对或批量改正。文章来源：https://www.toymoban.com/news/detail-423776.html
"""
GACOS大气改正python实现
"""
from make_correction import make_correction
import os
import numpy as np
# ---------------------------------------------------单干涉对改正---------------------------------------------------------
# dir = r'' # 绝对路径
# ifg_name = r'' # 干涉对文件名称
# ztd1_name = r''
# ztd2_name = r''
# ref_lat = 0
# ref_lon = 0
# wavelength = 0.055165
# unit_phase = 'P'
# unit_out = 'P'
# fln_elev = r''
# elev = 90
# isplanar = 1
# iswrap = 0
# make_correction(dir=dir, ifg_name=ifg_name, ztd1_name=ztd1_name, ztd2_name=ztd2_name, ref_lat=ref_lat, ref_lon=ref_lon,
# wavelength=wavelength, unit_phase=unit_phase, unit_out=unit_out, fln_elev=fln_elev, elev=elev,
# isplanar=isplanar, iswrap=iswrap)
# ----------------------------------------------------------批量改正-----------------------------------------------------
# ---------------------------------------需要保证文件夹下的数据日期是按照时间升序排列-------------------------------------------
dir = r'' # 绝对路径或是相对路径
save_dir = r'' # 保存文件的路径
process_dir = r'' # 中间处理数据结果的路径
ref_lat = 0
ref_lon = 0
wavelength = 0.055165 # phase wavelength,S1=0.055165
unit_phase = 'p' # interferogram unit,choose from p,m,c for phase, meter, centimeter
unit_out = 'p' # output unit,choose from p,m,c for phase, meter, centimeter
fln_elev = r'' # elevation angle file
elev = 90 # if elevation angle file not found, use this value instead,unit in degree
isplanar = 1 # if remove planar trend, 1 for yes, 0 for no
iswrap = 0 # if rewrap when plotting, note that the result file will not be wrapped
# 获取差分相位数据的个数---根据文件是否有后缀名来统计
# 获取文件列表
list_dir = os.listdir(dir)
# 获取差分文件名
diff = []
for i in range(len(list_dir)):
if os.path.splitext(list_dir[i])[-1] == '':
diff.append(list_dir[i])
else:
pass
# 获取各个时间的ztd名
single_diff = []
for i in range(len(diff)):
data = diff[i].split('-', 1)
single_diff.append(data)
# 去除重复名
ztd_name_list = np.unique(single_diff)
for i in range(len(diff)):
ztd1, ztd2 = diff[i].split('-', 1)
index1 = ztd_name_list.tolist().index(ztd1)
index2 = ztd_name_list.tolist().index(ztd2)
print(index1, index2)
ifg_file = ztd_name_list[index1] + '-' + ztd_name_list[index2]
dd = ztd_name_list[index1] + '.ztd'
make_correction(dir=dir, procecss_dir=process_dir, save_dir=save_dir, ifg_name=ifg_file, ztd1_name=ztd_name_list[index1] + '.ztd',
ztd2_name=ztd_name_list[index2] + '.ztd', ref_lat=ref_lat, ref_lon=ref_lon, wavelength=wavelength,
unit_phase=unit_phase, unit_out=unit_out, fln_elev=fln_elev, elev=elev, isplanar=isplanar,
iswrap=iswrap)
print()
8.ztd转los
由于gacos数据是基于ztd的，因此需要将其转到los向
注：这里入射角数据是对应研究区的矩阵形式，是一个入射角矩阵，由于sarscape没有生成这个文件，我这里简化了（直接将卫星入射角设置为矩阵），若能够得到入射角矩阵，下面的代码更改为element-wise形式即可文章来源地址https://www.toymoban.com/news/detail-423776.html
import struct
import numpy as np
import matplotlib.pyplot as plt
import os
def xshow(filename, nx, nz):
f = open(filename, "rb")
pic = np.zeros((nx, nz))
for i in range(nx):
for j in range(nz):
data = f.read(4)
elem = struct.unpack("f", data)[0]
pic[i][j] = elem
# print(pic)
f.close()
return pic
def ztd_los(dir, cols, lines, incidence, out_dir):
file = dir
cols = cols
lines = lines
file_list = os.listdir(file)
select_file_list = []
incidence = incidence
for i in range(len(file_list)):
select_file_list.append(file_list[i][0:17])
# 去除重复值
# updata_data_list = list(set(select_file_list))
updata_data_list = []
for i in range(len(select_file_list)):
if select_file_list[i] not in updata_data_list:
updata_data_list.append(select_file_list[i])
else:
pass
for i in range(len(updata_data_list)):
# set out path
dztd_los_path = out_dir + '\\' + updata_data_list[i] + '-dztd-los'
corrected_los_path = out_dir + '\\' + updata_data_list[i] + '-corrected-los'
# read data
dztd = xshow(filename=dir + '\\' + updata_data_list[i] + '-dztd', nx=cols, nz=lines)
corrected = xshow(filename=dir + '\\' + updata_data_list[i] + '-phase-dztd', nx=cols, nz=lines)
# flatten
dztd_flatten = dztd.flatten()
corrected_flatten = corrected.flatten()
# ztd-los
dztd_los = dztd_flatten / np.cos(incidence)
corrected_los = corrected_flatten / np.cos(incidence)
# 保存
with open(dztd_los_path, 'wb') as f:
for j in dztd_los:
a = struct.pack('f', j)
f.write(a)
print("dztd-los file Done")
print("-----------------------")
with open(corrected_los_path, 'wb') as f:
for j in corrected_los:
a = struct.pack('f', j)
f.write(a)
print("corrected-los file Done")
print("-----------------------")
print('The {} th GACOS result has been converted to LOS'.format(i+1))
print('============================================================')
# 可视化
dztd_los_data = xshow(dztd_los_path, cols, lines)
corrected_los_data = xshow(corrected_los_path, cols, lines)
# read geo file to get label
fig, ax = plt.subplots(1, 2)
ax1 = ax[0].imshow(dztd_los_data, cmap='jet')
cb1 = plt.colorbar(ax1, ax=ax[0], fraction=0.05)
ax2 = ax[1].imshow(corrected_los_data, cmap='jet')
cb2 = plt.colorbar(ax2, ax=ax[1], fraction=0.05)
cb1.ax.set_title('(rad)')
cb2.ax.set_title('(rad)')
fig.suptitle('GACOS results to LOS', fontsize=20, y=0.85)
ax[0].set_title('dztd-los')
ax[1].set_title('corrected-los')
plt.show()
到了这里，关于GACOS大气改正python实现的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！
