GACOS大气改正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]