sanya-dbms/processing/tools/Raster/MergeMetadata.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
)

# TODO 栅格数据文件夹
folder_path = 'E:\\projects\\遥感技术部需求\\批量生成元数据\\make\\data'
# TODO 元数据输出文件夹
output_path = 'E:\\projects\\遥感技术部需求\\批量生成元数据\\make\\metadata'
# TODO 元数据模板路径
templatepath = 'E:\\projects\\遥感技术部需求\\批量生成元数据\\make\\template\\template.xlsx'
# 获取元数据模板文件所在的目录
templatedirectory = os.path.dirname(os.path.abspath(templatepath))
# TODO 空间分析图层位置 默认在模板文件夹内有一个overlap.shp  点线面要素无所谓
overlapVectorPath = f"{templatedirectory}\\overlap.shp"
# 绑定需要从shp数据提取的信息
# data_gain_date 影像/数据获取时间
# data_resolution 影像分辨率/点云密度
overlapParams = {
    "data_resolution": "resolution",
    "data_gain_date": "gaindate"
}
# 通过文件路径来获取文件名和文件后缀
def getImageNumberAndFormat(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(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(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(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

# 遍历文件夹中的所有文件
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 = getImageNumberAndFormat(filepath)
        file_size = getFileSize(filepath)
        if dataset:
            # TODO 获取栅格数据的行数、列数
            grid_row = dataset.RasterYSize
            grid_col = dataset.RasterXSize
            # TODO 格网单元尺寸
            grid_size = 2
            # 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 = getImageLeftTop(dataset)
            # 获取地理信息
            geo_transform = dataset.GetGeoTransform()
            # 获取投影信息
            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 = spatial_ref.GetProjParm(osr.SRS_PP_CENTRAL_MERIDIAN)
                # TODO 椭球长半径
                ellipsoid_radius = 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")

        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": cornerpoint["lon"],
            "cornerpoint_lat": 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
        }
        # TODO 空间叠加
        for key in overlapParams:
            value = overlapParams[key]
            largest_feature = get_largest_overlapping_feature(raster_path=filepath, shp_path=overlapVectorPath)
            field_value = largest_feature[value]
            metadatadict[key] = field_value
        # 利用模板生成excel
        workbook = openpyxl.load_workbook(templatepath)
        sheet = workbook.active
        # 定义替换数据的字典
        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:
                            # 替换占位符
                            cell.value = cell.value.replace(placeholder, str(replacement))
        # 保存修改后的文件
        workbook.save(f'{output_path}\\{image_number}.xlsx')