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

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

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

from osgeo import gdal, osr
import openpyxl
from utils import *

# 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"
}
# 遍历文件夹中的所有文件
for filename in os.listdir(folder_path):
    filepath = os.path.join(folder_path, filename)
    # 检查文件扩展名
    if filename.endswith('.img') or filename.endswith('.tif'):
        # filepath = 'E:\\projects\\遥感技术部需求\\批量生成元数据\\数据\\分幅成果\\DEM成果\\DEM\\NH48G053055DEMK.img'
        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 = processingNumericalValues(spatial_ref.GetProjParm(osr.SRS_PP_CENTRAL_MERIDIAN))
                # TODO 椭球长半径
                ellipsoid_radius = 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")

        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')