148 lines
5.4 KiB
Python
148 lines
5.4 KiB
Python
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)
|
||
|