#!/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')