sanya-dbms/processing/tools/Raster/makeMetadate.py

# -*- coding: utf-8 -*-

__author__ = 'wanger'
__description__ = '批量生成栅格数据元数据表格 EXCEL'
__date__ = '2024-11-25'
__copyright__ = '(C) 2024 by siwei'
__revision__ = '1.0'

import os
from osgeo import gdal, osr
import openpyxl
from qgis.core import (
    QgsRasterLayer,
    QgsVectorLayer,
    QgsSpatialIndex,
    QgsFeatureRequest,
    QgsGeometry
)
from qgis.PyQt.QtCore import NULL
from qgis.PyQt.QtCore import QCoreApplication
from qgis.core import (QgsProcessing,
                       QgsFeatureSink,
                       QgsExpression,
                       QgsFeatureRequest,
                       QgsVectorLayer,
                       QgsProcessingException,
                       QgsProcessingAlgorithm,
                       QgsProcessingParameterFeatureSource,
                       QgsProcessingParameterFolderDestination,
                       QgsProcessingParameterString,
                       QgsProcessingParameterFile,
                       QgsProcessingParameterFeatureSink)
from qgis import processing


class MakeMetadataProcessingAlgorithm(QgsProcessingAlgorithm):
    rasterfolder = 'rasterfolder'
    template = 'template'
    outputfolder = 'outputfolder'
    # 图幅xy标识位取值 5:8,8:11
    splittingnum = 'splittingnum'
    splittingnumdefault = '5:8,8:11'

    def tr(self, string):
        return QCoreApplication.translate('Processing', string)

    def createInstance(self):
        return MakeMetadataProcessingAlgorithm()

    def name(self):
        return 'makeMetadata'

    def displayName(self):
        return self.tr('生成栅格元数据')

    def group(self):
        return self.tr('栅格数据')

    def groupId(self):
        return 'metadata'

    def shortHelpString(self):
        return self.tr("将选择文件夹内的所有栅格数据根据元数据模板动态生成元数据表。")

    def initAlgorithm(self, config=None):
        self.addParameter(QgsProcessingParameterFolderDestination(self.rasterfolder, '分幅数据文件夹'))
        self.addParameter(QgsProcessingParameterFile(self.template, '模板文件', extension='xlsx'))
        self.addParameter(QgsProcessingParameterFolderDestination(self.outputfolder, '元数据输出文件夹'))
        self.addParameter(
            QgsProcessingParameterString(self.splittingnum, self.tr('图幅xy标识位'), self.splittingnumdefault))

    # 通过文件路径来获取文件名和文件后缀
    def getImageNumberAndFormat(self, filepath):
        # 获取文件名（包括后缀）
        file_name = os.path.basename(filepath)
        # 获取文件名（不包括后缀）
        file_name_no_ext = os.path.splitext(file_name)[0]
        # 获取文件后缀
        file_ext = os.path.splitext(file_name)[1][1:]
        return file_name_no_ext, file_ext

    # 通过文件路径获取文件大小
    def getFileSize(self, filepath):
        try:
            file_size = os.path.getsize(filepath)
            size_mb = file_size / (1024 * 1024)
            return f"{size_mb:.2f}"
        except FileNotFoundError:
            print("文件未找到，请检查路径是否正确！")
            return 0

    # 获取左上角点
    def getImageLeftTop(self, dataset):
        # 获取地理变换参数
        geotransform = dataset.GetGeoTransform()
        if not geotransform:
            raise ValueError("无法获取地理变换参数！")

        # 左上角像素的坐标（投影坐标系）
        x_min = geotransform[0]
        y_max = geotransform[3]

        # 获取栅格的空间参考（投影坐标系）
        spatial_ref = osr.SpatialReference()
        spatial_ref.ImportFromWkt(dataset.GetProjection())
        # 转换为地理坐标系（经纬度）
        if spatial_ref.IsProjected():
            # spatial_ref_geo = spatial_ref.CloneGeogCS()
            wgs84_crs = osr.SpatialReference()
            wgs84_crs.ImportFromEPSG(4326)
            transform = osr.CoordinateTransformation(spatial_ref, wgs84_crs)
            lon, lat, _ = transform.TransformPoint(x_min, y_max)
        else:
            lon, lat = x_min, y_max  # 如果已经是地理坐标系，直接使用
        return {
            "lon": lat,
            "lat": lon,
            "y": int(x_min),
            "x": int(y_max)
        }

    # 通过栅格数据与矢量数据叠加分析  获取相交面积最大的一块  返回要素
    def get_largest_overlapping_feature(self, raster_path, shp_path):
        """
        查询与栅格相交的矢量要素中面积最大的一块，并返回指定字段的值。

        :param raster_path: 栅格文件路径
        :param shp_path: 矢量文件路径
        :return: 面积最大的要素
        """
        # 加载栅格图层
        raster_layer = QgsRasterLayer(raster_path, "Raster Layer")
        if not raster_layer.isValid():
            raise Exception(f"栅格文件无效: {raster_path}")
        # 加载矢量图层
        vector_layer = QgsVectorLayer(shp_path, "Vector Layer", "ogr")
        if not vector_layer.isValid():
            raise Exception(f"SHP文件无效: {shp_path}")
        # 获取栅格的范围
        raster_extent = raster_layer.extent()
        # 创建空间索引以加速查询
        spatial_index = QgsSpatialIndex(vector_layer.getFeatures())
        # 查询与栅格范围相交的要素ID
        intersecting_ids = spatial_index.intersects(raster_extent)
        # 遍历这些要素并计算相交区域的面积
        max_area = 0
        largest_feature = None
        for fid in intersecting_ids:
            feature = vector_layer.getFeature(fid)
            intersection = feature.geometry().intersection(QgsGeometry.fromRect(raster_extent))
            # 检查是否有有效相交区域
            if intersection.isEmpty():
                continue
            # 计算相交区域的面积
            area = intersection.area()
            # 找到面积最大的要素
            if area > max_area:
                max_area = area
                largest_feature = feature
        if largest_feature:
            return largest_feature
        else:
            return None

    # 替换字符串  type (first,last)
    def replacestr(self, string, oldvalue, newvalue, type):
        oldvalue = str(oldvalue)
        newvalue = str(newvalue)
        if type == "first":
            return string.replace(oldvalue, newvalue, 1)
        elif type == "last":
            # 查找最后一次出现的 oldvalue 的位置
            index = string.rfind(oldvalue)
            if index != -1:
                # 用切片将 last occurrence 替换
                return string[:index] + newvalue + string[index + len(oldvalue):]

    # 获取接边情况
    def getSplitting(self, image_number, image_format, imagex, imagey, length, folder_path, selectfeature):
        east = f"{self.replacestr(string=image_number, oldvalue=str(imagey).zfill(length), newvalue=str(imagey + 1).zfill(length), type='last')}.{image_format}"
        west = f"{self.replacestr(string=image_number, oldvalue=str(imagey).zfill(length), newvalue=str(imagey - 1).zfill(length), type='last')}.{image_format}"
        north = f"{self.replacestr(string=image_number, oldvalue=str(imagex).zfill(length), newvalue=str(imagex - 1).zfill(length), type='first')}.{image_format}"
        south = f"{self.replacestr(string=image_number, oldvalue=str(imagex).zfill(length), newvalue=str(imagex + 1).zfill(length), type='first')}.{image_format}"
        print(selectfeature["east"])
        print(selectfeature["west"])
        print(selectfeature["north"])
        print(selectfeature["south"])
        return {
            "east": selectfeature["east"] if selectfeature["east"] != NULL else self.hasFile(east, folder_path),
            "west": selectfeature["west"] if selectfeature["west"] != NULL else self.hasFile(west, folder_path),
            "north": selectfeature["north"] if selectfeature["north"] != NULL else self.hasFile(north, folder_path),
            "south": selectfeature["south"] if selectfeature["south"] != NULL else self.hasFile(south, folder_path),
        }

    # 判断有没有文件
    def hasFile(self, filepath, folder_path):
        if os.path.exists(f"{folder_path}\\{filepath}"):
            return "已接"
        else:
            return "未接"

    # 处理数值
    def processingNumericalValues(self, value):
        print(value)
        # return value
        if isinstance(value, str):
            try:
                value = float(value)
            except ValueError:
                print(f"无法将 '{value}' 转换为浮动数")
                return None
        if value.is_integer():
            return int(value)
        else:
            return round(value, 2)

    # 经纬度转度分秒
    def decimal_to_dms(self, deg):
        # 计算度（Degrees）
        degrees = int(deg)
        # 计算分（Minutes），并截取小数部分
        minutes = int((deg - degrees) * 60)
        # 计算秒（Seconds），保留两位小数
        seconds = round((deg - degrees - minutes / 60) * 3600, 2)
        return f"{degrees}{minutes}{seconds}"

    # 运行
    def processAlgorithm(self, parameters, context, feedback):
        # 获取模板文件夹
        template_file = self.parameterAsString(parameters, self.template, context)
        # 栅格数据文件夹
        raster_folder = self.parameterAsString(parameters, self.rasterfolder, context)
        # 元数据生成目录
        output_folder = self.parameterAsString(parameters, self.outputfolder, context)
        # 通过图号取接边数据标识位
        splittingnum = self.parameterAsString(parameters, self.splittingnum, context)
        print(template_file, raster_folder, output_folder)
        # TODO 栅格数据文件夹
        folder_path = raster_folder
        # TODO 元数据输出文件夹
        output_path = output_folder
        # TODO 元数据模板路径
        templatepath = template_file
        # 获取元数据模板文件所在的目录
        templatedirectory = os.path.dirname(os.path.abspath(templatepath))
        # TODO 空间分析图层位置 默认在模板文件夹内有一个overlap.shp  点线面要素无所谓
        overlapVectorPath = f"{templatedirectory}\\overlap.shp"
        # TODO 网格数据查询 默认在模板文件夹内有一个fishnet.shp  点线面要素无所谓
        fishVectorPath = f"{templatedirectory}\\fishnet.shp"
        fishlayer = QgsVectorLayer(fishVectorPath, "Fish Layer", "ogr")
        fish_field_name = "imagenum"
        # 绑定需要从shp数据提取的信息
        # data_gain_date 影像/数据获取时间
        # data_resolution 影像分辨率/点云密度
        overlapParams = {
            "data_resolution": "resolution",
            "data_gain_date": "gaindate",
            "source_type": "sourcetype"
        }
        # 记录总计生成了多少个元数据表
        totalNumber = 0
        # 遍历文件夹中的所有文件
        for filename in os.listdir(folder_path):
            filepath = os.path.join(folder_path, filename)
            # 检查文件扩展名
            if filename.endswith('.img') or filename.endswith('.tif'):
                dataset = gdal.Open(filepath)
                # TODO 读取图号
                image_number, image_format = self.getImageNumberAndFormat(filepath)
                file_size = self.getFileSize(filepath)
                # 通过图号获取网格要素
                sql_query = f"strpos('{image_number}', \"{fish_field_name}\") > 0"
                expression = QgsExpression(sql_query)
                request = QgsFeatureRequest(expression)
                # 获取查询结果
                matching_features = [f for f in fishlayer.getFeatures(request)]
                selectfeature = matching_features[0]
                if dataset:
                    # TODO 获取栅格数据的行数、列数
                    grid_row = dataset.RasterYSize
                    grid_col = dataset.RasterXSize
                    # TODO波段数
                    band_size = dataset.RasterCount
                    # TODO 无效值 list 用逗号隔开
                    nodata_list = []
                    # TODO 计算高程记录的小数点位数
                    decimal_number = None
                    for band_index in range(1, dataset.RasterCount + 1):
                        band = dataset.GetRasterBand(band_index)
                        nodata_list.append(band.GetNoDataValue())
                        stats = band.GetStatistics(True, True)  # 获取统计信息
                        min_value, max_value = stats[0], stats[1]
                        if min_value is not None and max_value is not None:
                            # 获取高程小数位数的最大可能值
                            decimal_number = max(len(str(min_value).split(".")[-1]),
                                                 len(str(max_value).split(".")[-1])) if '.' in str(
                                min_value) else 0
                    # TODO 获取栅格数据左上角标
                    cornerpoint = self.getImageLeftTop(dataset)
                    # 获取地理信息
                    geo_transform = dataset.GetGeoTransform()
                    # TODO 格网单元尺寸
                    x_resolution = abs(geo_transform[1])  # 像素宽度
                    y_resolution = abs(geo_transform[5])  # 像素高度
                    grid_size = f"{self.processingNumericalValues(x_resolution)}*{self.processingNumericalValues(y_resolution)}"
                    # grid_size = f"{x_resolution}*{y_resolution}"
                    # 获取投影信息
                    projection = dataset.GetProjection()
                    # 获取栅格数据的投影信息
                    # 创建 SpatialReference 对象并加载投影信息
                    spatial_ref = osr.SpatialReference()
                    spatial_ref.ImportFromWkt(projection)
                    # TODO 投影参数集合
                    central_meridian = None
                    ellipsoid_radius = None
                    ellipsoid_flat = None
                    geodetic_datum = None
                    projection_name = None
                    divided_into_zones = "无分带"
                    with_number = "无代号"
                    plane_coordinate_unit = "米"
                    altitude_name = "正常高"
                    altitude_benchmark = "1985国家高程基准"
                    if spatial_ref.IsProjected():
                        # TODO 中央子午线
                        central_meridian = self.processingNumericalValues(
                            spatial_ref.GetProjParm(osr.SRS_PP_CENTRAL_MERIDIAN))
                        # TODO 椭球长半径
                        ellipsoid_radius = self.processingNumericalValues(spatial_ref.GetSemiMajor())  # 长半轴
                        # TODO 椭球扁率
                        ellipsoid_flat = f"1/{spatial_ref.GetInvFlattening():.3f}"  # 反扁率（1/扁率）
                        # TODO 所采用大地基准
                        geodetic_datum = spatial_ref.GetAttrValue("DATUM")
                        # TODO 地图投影名称
                        projection_name = spatial_ref.GetAttrValue("PROJECTION")
                        # TODO 分带方式（适用于 UTM）
                        if geodetic_datum == "China_2000":
                            if projection.__contains__("-degree"):
                                number = projection.split("-degree")[0][-1:]
                                divided_into_zones = f'{number}度带'
                                with_number = int(central_meridian / int(number))
                        elif spatial_ref.IsProjected() and "UTM" in projection:
                            zone = spatial_ref.GetUTMZone()
                            divided_into_zones = "UTM分带" if zone else "其他"
                            with_number = zone if zone else "N/A"
                        # TODO 平面坐标单位
                        plane_coordinate_unit = spatial_ref.GetLinearUnitsName()
                        linear_units = spatial_ref.GetLinearUnits()  # 单位与米的换算关系
                        # 高程系统（如果有）
                        altitude_name = spatial_ref.GetAttrValue("VERT_CS")
                        altitude_benchmark = spatial_ref.GetAttrValue("VERT_DATUM")
                # 获取四向接边情况
                numlist = splittingnum.split(",")
                length = int(numlist[0].split(":")[1]) - int(numlist[0].split(":")[0])
                imagex = int(image_number[int(numlist[0].split(":")[0]):int(numlist[0].split(":")[1])])
                imagey = int(image_number[int(numlist[1].split(":")[0]):int(numlist[1].split(":")[1])])
                splitting = self.getSplitting(image_number, image_format, imagex, imagey, length, folder_path,
                                              selectfeature)
                metadatadict = {
                    "image_number": image_number,
                    "data_format": image_format,
                    "data_volume": file_size,
                    "grid_size": grid_size,
                    "grid_row": grid_row,
                    "grid_col": grid_col,
                    "decimal_number": decimal_number,
                    "nodata_number": ",".join(map(str, nodata_list)),
                    "cornerpoint_lon": self.decimal_to_dms(cornerpoint["lon"]),
                    "cornerpoint_lat": self.decimal_to_dms(cornerpoint["lat"]),
                    "cornerpoint_x": cornerpoint["x"],
                    "cornerpoint_y": cornerpoint["y"],
                    "ellipsoid_radius": ellipsoid_radius,
                    "ellipsoid_flat": ellipsoid_flat,
                    "geodetic_datum": "2000国家大地坐标系" if geodetic_datum == "China_2000" else geodetic_datum,
                    "projection_name": "高斯-克吕格投影" if projection_name == "Transverse_Mercator" else projection_name,
                    "central_meridian": central_meridian,
                    "divided_into_zones": divided_into_zones,
                    "with_number": with_number,
                    "plane_coordinate_unit": "米" if plane_coordinate_unit == "metre" else plane_coordinate_unit,
                    "altitude_name": altitude_name,
                    "altitude_benchmark": altitude_benchmark,
                    "join_edge_w": splitting["west"],
                    "join_edge_e": splitting["east"],
                    "join_edge_n": splitting["north"],
                    "join_edge_s": splitting["south"],
                    "land_forms": selectfeature["landforms"]
                }
                for field_name in selectfeature.fields().names():
                    metadatadict[field_name] = selectfeature[field_name]
                print("栅格数据参数获取完成")
                # TODO 空间叠加
                for key in overlapParams:
                    value = overlapParams[key]
                    largest_feature = self.get_largest_overlapping_feature(raster_path=filepath,
                                                                           shp_path=overlapVectorPath)
                    field_value = largest_feature[value]
                    metadatadict[key] = field_value
                print("开始生成报告")
                # print(metadatadict)
                # 利用模板生成excel
                workbook = openpyxl.load_workbook(templatepath)
                print(f"{templatepath}打开成功")
                sheet = workbook.active
                print(f"sheet打开成功")
                # 定义替换数据的字典
                replacements = {}
                for attr in metadatadict:
                    key = f"<{attr}>"
                    value = metadatadict[attr]
                    replacements[key] = value
                # 遍历单元格，替换占位符
                for row in sheet.iter_rows():
                    for cell in row:
                        if isinstance(cell.value, str):  # 确保单元格内容是字符串
                            for placeholder, replacement in replacements.items():
                                if placeholder in cell.value:
                                    # 替换占位符
                                    try:
                                        cell.value = cell.value.replace(placeholder, str(replacement))
                                    except SomeException as e:
                                        print(e)
                # 保存修改后的文件
                print(output_path)
                workbook.save(f'{output_path}\\{image_number}.xlsx')
                totalNumber = totalNumber + 1
                print(f"元数据生成完成:{image_number}")
        print(f"共生成{totalNumber}个元数据表。")
        return {
            "状态": "成功",
            "任务": f"共生成{totalNumber}个元数据表。"
        }