#!/usr/bin/env
|
# -*- coding: utf-8 -*-
|
|
|
import os
|
import sys
|
import time
|
|
sys.path.append(r"C:\Program Files\QGIS 3.16\apps\qgis-ltr\python\plugins")
|
import math
|
import argparse
|
from qgis.core import *
|
from qgis.gui import *
|
from qgis.PyQt.QtGui import *
|
from qgis.PyQt.QtCore import *
|
from qgis.PyQt.QtWidgets import *
|
from processing.core.Processing import Processing
|
|
|
# 获取完整路径
|
def get_full_path():
|
return os.path.split(sys.argv[0])[0]
|
|
|
# 获取参数
|
def get_args():
|
print("argv = ", sys.argv[1:])
|
parser = argparse.ArgumentParser(description='ArgUtils')
|
parser.add_argument("-src", type=str, default=get_full_path(), required=False)
|
parser.add_argument("-qgz", type=str, default=r"xyz.qgz", required=False)
|
parser.add_argument("-file", type=str, default=r"D:\xyz\dem\dem.txt", required=False)
|
parser.add_argument("-out", type=str, default=r"D:\xyz\dem\out.tif", required=False)
|
|
return parser.parse_args()
|
|
|
# 读取文本文件
|
def read_txt(filePath):
|
f = open(filePath, encoding="utf-8")
|
lines = f.readlines()
|
f.close()
|
|
for i in range(0, len(lines)):
|
lines[i] = lines[i].replace('\n', '')
|
|
return lines
|
|
|
# 加载图层
|
def load_layers(prj, args):
|
for layer in prj.mapLayers().values():
|
prj.removeMapLayer(layer)
|
|
args.authid = None
|
lines = read_txt(args.file)
|
for i in range(0, len(lines)):
|
if len(lines[i]) == 0:
|
continue
|
|
print("layer_" + str(i) + ": " + lines[i])
|
layer = QgsRasterLayer(lines[i], "layer_" + str(i))
|
if not layer.isValid() or layer.crs() is None:
|
print("layer_" + str(i) + ": failed to load!")
|
continue
|
|
if len(prj.mapLayers()) == 0:
|
args.authid = layer.crs().authid()
|
print("authid: " + args.authid)
|
prj.setCrs(layer.crs())
|
|
# for j in range(1, layer.bandCount() + 1):
|
# try:
|
# layer.dataProvider().setNoDataValue(j, 0)
|
# except Exception as e:
|
# print(e)
|
|
prj.addMapLayer(layer)
|
|
|
# 获取参数
|
def get_ops(prj, args):
|
layers = []
|
for layer in prj.mapLayers().values():
|
layers.append(layer.name())
|
|
ops = {
|
'INPUT': layers,
|
'DATA_TYPE': 5,
|
'OUTPUT': args.out,
|
# 'NODATA_INPUT': 0,
|
# 'NODATA_OUTPUT': 0,
|
'OPTIONS': ''
|
}
|
print(ops)
|
|
return ops
|
|
|
# 合并
|
def merge(prj, args):
|
import processing
|
|
ops = get_ops(prj, args)
|
processing.run("gdal:merge", ops)
|
|
|
# 获取中心点坐标
|
def get_coord(prj, args):
|
layer = QgsRasterLayer(args.out, "layer_tif")
|
rect = layer.extent()
|
x = (rect.xMinimum() + rect.xMaximum()) / 2
|
y = (rect.yMinimum() + rect.yMaximum()) / 2
|
|
sid = layer.crs().authid()
|
if sid != "EPSG:4490" and sid != "EPSG:4326":
|
sid = "EPSG:4490" if layer.crs().toWkt().find("CGCS2000") > -1 else "EPSG:4326"
|
transform = QgsCoordinateTransform(QgsCoordinateReferenceSystem(layer.crs().authid()), QgsCoordinateReferenceSystem(sid), prj)
|
|
p = QgsPoint(x, y)
|
p.transform(transform)
|
x = p.x()
|
y = p.y()
|
|
coord = str(x) + ", " + str(y) + ", " + str(sid)
|
print(coord)
|
|
f = open(args.out.replace(".tif", "_cs.txt"), 'w')
|
f.write(coord)
|
f.close()
|
|
|
# 初始化
|
def init():
|
# QgsApplication.setPrefixPath("C:\Program Files\QGIS 3.16", True)
|
qgs = QgsApplication([], False)
|
qgs.initQgis()
|
Processing.initialize()
|
|
args = get_args()
|
prj = QgsProject.instance()
|
prj.read(os.path.join(args.src, args.qgz))
|
# prj.read(args.qgz)
|
print("FileName: " + prj.fileName())
|
|
load_layers(prj, args)
|
merge(prj, args)
|
get_coord(prj, args)
|
|
qgs.exitQgis()
|
|
|
# main
|
if __name__ == '__main__':
|
timer = time.time()
|
init()
|
|
print(f'耗时:{time.time() - timer:.2f}s')
|
