Python Bouguet校正算法 python大气校正

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.主程序

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

"""
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形式即可

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()