microproduct-l-sar/atmosphericDelay-L-SAR/read_post.py

import glob
import os
import pandas as pd
from pathlib import Path
import csv
import numpy as np
from osgeo import gdal
from xml.etree.ElementTree import ElementTree, Element
import shutil


def searchfile(source, pattern='*', recursive=False):

    if not isinstance(source, Path):
        source = Path(source)
    if source.is_dir() and source.exists():
        if recursive:
            res = list(source.rglob(pattern))
        else:
            res = list(source.glob(pattern))
    else:
        print(f"{source} is a file path or folder path is not exists!")
        return []   # False
    file_list = [i for i in res if i.is_file()]
    return file_list


def read_posTrp(file_dir):
    psoT = []
    lon = []
    lat = []
    file_name = os.path.join(file_dir, "posTrp1.csv")
    posTrp_file = searchfile(file_dir, '*.posTrp', recursive=True)
    for file in posTrp_file:
        pos = open(file, 'r')
        lines = pos.readlines()
        for line in lines:
            a = line.split()
            if len(a) == 11:
                if float(a[3]) <= 0 :
                    continue
                else:
                    lon = a[4][0:14]
                    x = [a[0], a[1] + ' ' + a[2], a[3], lon, a[8], a[9], a[10]]
                # if float(a[8]) == 0:
                #     continue
                # else:
                #     x = [a[0], a[1] + ' ' + a[2], a[3], a[4], a[8], a[9], a[10]]
            elif len(a) == 12:
                if float(a[3]) <= 0 or float(a[4]) <= 0:
                    continue
                else:
                    x = [a[0], a[1] + ' ' + a[2], a[3], a[4], a[9], a[10], a[11]]
                # if float(a[9]) == 0:
                #     continue
                # else:
                #     x = [a[0], a[1] + ' ' + a[2], a[3], a[4], a[9], a[10], a[11]]
            else:
                if float(a[3]) <= 0 or float(a[4]) <= 0:
                    continue
                else:
                    x = [a[0], a[1] + ' ' + a[2], a[3], a[4], a[6], a[7], a[8]]

                lon.append(float(a[4]))
                lat.append(float(a[3]))
            psoT.append(x)
    with open(file_name, 'w', newline='') as writ:
        writer = csv.writer(writ)
        writer.writerow(['station', 'time', 'lat', 'lon', 'delay_all', 'delay_dry', 'delay_wet'])
        writer.writerows(psoT)

    # print(np.max(lon))
    # print(np.min(lon))
    # print(np.max(lat))
    # print(np.min(lat))


def create_para(dem_file):
    para_dict = {}
    gdal.AllRegister()
    basename = os.path.basename(dem_file)
    names = basename.split("_")
    name = basename.split(".")
    para_dict.update({"name": name[0]})
    para_dict.update({"datatype": names[0]})
    dataset = gdal.Open(str(dem_file))
    trans = dataset.GetGeoTransform()
    # [(左上角像元的位置), (像元的宽度), (如果影像是指北的，0),(左上角像元的位置),（如果影像是指北的，0）, （像元的高度）]
    weight = dataset.RasterXSize
    height = dataset.RasterYSize

    TopLeftLatitude = trans[3]
    TopLeftLongitude = trans[0]
    TopRightLatitude = trans[3]
    TopRightLongitude = trans[0] + trans[1] * weight
    BottomleftLatitude = trans[3] + trans[5] * height
    BottomleftLongitude = trans[0]
    BottomrightLatitude = trans[3] + trans[5] * height
    BottomrightLongitude = trans[0] + trans[1] * weight
    para_dict.update({"imageinfo_TopLeftLatitude": TopLeftLatitude})
    para_dict.update({"imageinfo_TopLeftLongitude": TopLeftLongitude})
    para_dict.update({"imageinfo_TopRightLatitude": TopRightLatitude})
    para_dict.update({"imageinfo_TopRightLongitude": TopRightLongitude})
    para_dict.update({"imageinfo_BottomleftLatitude": BottomleftLatitude})
    para_dict.update({"imageinfo_BottomleftLongitude": BottomleftLongitude})
    para_dict.update({"imageinfo_BottomrightLatitude": BottomrightLatitude})
    para_dict.update({"imageinfo_BottomrightLongitude": BottomrightLongitude})
    return para_dict


def create_xml(para_dict, model_xml, out_path):
    shutil.copy(model_xml, out_path)
    tree = ElementTree()
    tree.parse(model_xml)  # 影像头文件
    root = tree.getroot()
    productinfo = tree.getroot()
    for key, value in para_dict.items():
        if key.split("_")[0] != "imageinfo":
            productinfo.find(key).text = str(value)
        elif key.split("_")[0] == "imageinfo":
            imageinfo = productinfo.find("spatRepInfo")
            if key.split("_")[1] in ["TopLeftLatitude", "TopLeftLongitude", "TopRightLatitude", "TopRightLongitude",
                                     "BottomleftLatitude", "BottomleftLongitude", "BottomrightLatitude",
                                     "BottomrightLongitude"]:
                imageinfo.find(key.split("_")[1]).text = str(value)
                pass
    tree.write(out_path, encoding="utf-8", xml_declaration=True)



def create_demxml(dem_dir, model_xml):
    dem_file = searchfile(dem_dir, '*.tif', recursive=True)
    for file in dem_file:
        para_dict = create_para(file)
        out_path = str(file).replace(".tif",".xml")
        if os.path.exists(out_path):
            os.remove(out_path)
        create_xml(para_dict, model_xml, out_path)


if __name__ == '__main__':
    file_dir = r'F:\干涉大气延迟校正\GNSS\2022'
    read_posTrp(file_dir)
    # dem_dir = r'F:\MicroWorkspace\mico_datas\Micro\astergdem2'
    # model_xml = r'F:\MicroWorkspace\mico_datas\Micro\model.xml'
    # create_demxml(dem_dir, model_xml)