GACOS大气改正python实现

该代码共有8个部分。分别是数据读取、头文件读取、ztd数据裁剪、趋势项去除、相位包裹、make correction、ztd转los以及主程序

我为什么要写这个代码呢？原因有二：
首先，因为我个人不喜欢使用matlab，用着不顺手。
其次，我想锻炼一下自己的编程能力，所以就着手了并完成了这样一个问题。
注：代码有任何问题，请私信我，接受建议，谢谢！！！

注：转载请标明出处，谢谢！！！

1.数据读取

常规的二进制文件读取。可以自己先调试数据是否读取正确，可视化一下，看是否正常。代码命名为read_binary.py

import numpy as np
import structdef 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] = elemf.close()return pic

2.ztd头文件读取

读取的就是下载的gacos头文件数据（.rsc）。代码命名为read_rsc.py

def read_rsc(file):file_path = filewith 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:passreturn 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 = filewith 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:passreturn width, length, x_first, y_first, x_step, y_stepdef cut_image(dir, out_path, filename, maxlat, len_l, minlon, wid_l, outfilename):file_path = dir + '\\' + filenamersc_file_path = dir + '\\' + filename + '.rsc'cut_outpath = out_path + '\\' + filename + outfilenamewidth, length, x_first, y_first, x_step, y_step = read_rsc(rsc_file_path)# read ztd dataztd_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 * ilon[i, j] = x_first + x_step * jindex_maxlat = round((maxlat - y_first) / y_step)index_minlat = index_maxlat + len_l - 1index_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 pltdef remove_planar(corrected):num = 0sumx = 0sumy = 0sumxy = 0sumx2 = 0sumy2 = 0sumsx = 0sumsy = 0sums = 0# data = np.array(corrected)data = corrected.flatten()  # 为了获取index，所以将数据展平，后续操作还是使用corrected# data = correctedzero_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与0continuegridx = jgridy = iphase = corrected[i][j]sumx = sumx + gridxsumy = sumy + gridysumx2 = sumx2 + gridx * gridxsumy2 = sumy2 + gridy * gridysumxy = sumxy + gridx * gridysumsx = sumsx + phase * gridxsumsy = sumsy + phase * gridysums = sums + phasenum = num + 1A = np.zeros((3, 3))A[0, 0] = sumx2A[0, 1] = sumxyA[0, 2] = sumxA[1, 0] = sumxyA[1, 1] = sumy2A[1, 2] = sumyA[2, 0] = sumxA[2, 1] = sumyA[2, 2] = numL = np.zeros(3)L[0] = sumsxL[1] = sumsyL[2] = sumsL = 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:passgridx = jgridy = itdata[i][j] = corrected[i][j] - R[0]*gridx - R[1] * gridy -R[2]new_tdata = tdata.flatten()new_tdata[zero_index] = 0return 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 pltdef 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].indexsize = min_index.sizeoutput = dataoutput_flatten = output.flatten()while min_index.size != 0:# 数据均是数组--根据min_index的位置将其数值处理--不包括nanoutput_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.nanplt.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 operatordef 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 / 100else:raise ValueError('please check the read unit')# -------------------------------------------read ztd files from GACOS----------------------------------------------ztd1_data = cut_ztd1  # data type：ndarrayztd2_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 = 0ifg_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写入原始数据中，方便找到对应的indexsub_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'].valuesmin_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:passchange_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.0new_dztd_data = dztd_data * 100.0else: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.0dztd = 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'].valuesdztd_label = dztd_dataframe['dztd'].valuesfor i in range(len(phase_label)):if math.isnan(phase_label[i]):dztd_label[i] = np.nanwith 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.主程序

实现单干涉对或批量改正。

"""
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向

import struct
import numpy as np
import matplotlib.pyplot as plt
import osdef 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 picdef ztd_los(dir, cols, lines, incidence, out_dir):file = dircols = colslines = linesfile_list = os.listdir(file)select_file_list = []incidence = incidencefor 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:passfor i in range(len(updata_data_list)):# set out pathdztd_los_path = out_dir + '\\' + updata_data_list[i] + '-dztd-los'corrected_los_path = out_dir + '\\' + updata_data_list[i] + '-corrected-los'# read datadztd = 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)# flattendztd_flatten = dztd.flatten()corrected_flatten = corrected.flatten()# ztd-losdztd_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 labelfig, 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实现相关推荐

Github配置(git+vscode+python+jupyter)
①下载git 打开 git bash 工具的用户名和密码存储 $ git config --global user.name "Your Name" $ git config -- ...
【实验楼】python简明教程
①终端输入python进入欣赏完自己的杰作后,按 Ctrl + D 输入一个 EOF 字符来退出解释器,你也可以键入 exit() 来退出解释器. ②vim键盘快捷功能分布 ③这里需要注意如果程序中 ...
【Kaggle Learn】Python 5-8
五. Booleans and Conditionals Using booleans for branching logic x = True print(x) print(type(x))''' ...
【Kaggle Learn】Python 1-4
[Kaggle Learn]Python https://www.kaggle.com/learn/python 一. Hello, Python A quick introduction to Py ...
使用python愉快地做高数线代题目~
今天接触到了python,发现真是极易上手啊!对比c语言是什么鬼东西= = 诶,等下,看完教学文章发现TA在下面写了这句话如果做了前面的内容你可能已被吸引了,觉得c语言真的是废材! 不...不是的. ...
python 位运算与等号_Python 运算符
和大多数语言一样,Python也有很多运算符,并且运算符跟其他语言的运算符大同小异接下来一一介绍: 算术运算符: 运算符描述实例 +加 - 两个对象相加a+b的输出结果是30 -减 - 得到复数或者一 ...
python减小内存占用_如何将Python内存占用缩小20倍？
当程序执行过程中RAM中有大量对象处于活动状态时,可能会出现内存问题,特别是在对可用内存总量有限制的情况下. 下面概述了一些减小对象大小的方法,这些方法可以显著减少纯Python程序所需的RAM数量. ...
python中排序英文单词怎么写_Python实现对文件进行单词划分并去重排序操作示例...
本文实例讲述了Python实现对文件进行单词划分并去重排序操作.,具体如下: 文件名:test1.txt 文件内容: But soft what light through yonder window ...
python程序如何执行死刑图片_如何判断对象已死
已死的对象就是不可能被任何途径使用的对象,有以下几种方法判断一个对象是否已经死了: 引用计数给对象添加一个引用计数器,每当有一个地方引用他,计算器就加 1:当引用失效时,计数器减 1:任何时刻计数器 ...

GACOS大气改正python实现

GACOS大气改正python实现

该代码共有8个部分。分别是数据读取、头文件读取、ztd数据裁剪、趋势项去除、相位包裹、make correction、ztd转los以及主程序

注：转载请标明出处，谢谢！！！

1.数据读取

2.ztd头文件读取

3.ztd数据裁剪

4.趋势项去除

5.相位包裹

6.make correction

7.主程序

8.ztd转los

GACOS大气改正python实现相关推荐

最新文章

热门文章