Преглед изворни кода

栅格数据元数据生成

wanger пре 7 месеци
родитељ
комит
9e1cf48275

+ 115 - 0
processing/tools/Raster/FishNet.py

@@ -0,0 +1,115 @@
+# -*- coding: utf-8 -*-
+
+__author__ = 'wanger'
+__description__ = '处理渔网数据  判断网格是否为自由 增加标签'
+__date__ = '2024-11-25'
+__copyright__ = '(C) 2024 by siwei'
+__revision__ = '1.0'
+
+import os
+from qgis.PyQt.QtCore import QCoreApplication
+from qgis.core import (QgsProcessing,
+                       QgsVectorLayer,
+                       QgsFeatureSink,
+                       QgsProcessingException,
+                       QgsProcessingAlgorithm,
+                       QgsProcessingParameterFile,
+                       QgsProcessingParameterFeatureSource,
+                       QgsProcessingParameterFeatureSink)
+from qgis import processing
+
+
+class FishNetProcessingAlgorithm(QgsProcessingAlgorithm):
+    INPUT = 'INPUT'
+
+    def tr(self, string):
+        return QCoreApplication.translate('Processing', string)
+
+    def createInstance(self):
+        return FishNetProcessingAlgorithm()
+
+    def name(self):
+        return 'makeFishNet'
+
+    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(QgsProcessingParameterFile(
+            self.INPUT,
+            '网格数据',
+            extension='shp'  # 设置文件过滤器,只允许选择shp文件
+        ))
+
+    # 定义检查一个网格单元格上下左右是否有相邻的网格单元
+    def check_neighbors(self, feature, layer):
+        # 获取当前单元格的几何
+        geom = feature.geometry()
+        bbox = geom.boundingBox()  # 获取网格单元的包围盒(bounding box)
+        # 提取当前网格的中心点位置
+        center_x, center_y = geom.centroid().asPoint()
+        # 定义上下左右偏移量(单位为网格的大小,假设网格是正方形或矩形)
+        tolerance = 1e-5  # 设定一个微小的容忍误差,用于比较几何形状
+        neighbors = {
+            'north': None,
+            'south': None,
+            'west': None,
+            'east': None
+        }
+        for other_feature in layer.getFeatures():
+            if other_feature.id() == feature.id():
+                continue  # 跳过自己
+            other_geom = other_feature.geometry()
+            # 检查上、下、左、右的邻接
+            if abs(center_y - other_geom.centroid().asPoint().y()) < tolerance and other_geom.centroid().asPoint().x() < center_x:
+                neighbors['west'] = other_feature
+            elif abs(
+                    center_y - other_geom.centroid().asPoint().y()) < tolerance and other_geom.centroid().asPoint().x() > center_x:
+                neighbors['east'] = other_feature
+            elif abs(
+                    center_x - other_geom.centroid().asPoint().x()) < tolerance and other_geom.centroid().asPoint().y() < center_y:
+                neighbors['south'] = other_feature
+            elif abs(
+                    center_x - other_geom.centroid().asPoint().x()) < tolerance and other_geom.centroid().asPoint().y() > center_y:
+                neighbors['north'] = other_feature
+        return neighbors
+
+    # 执行
+    def processAlgorithm(self, parameters, context, feedback):
+        shp_file_path = self.parameterAsString(parameters, self.INPUT, context)
+        # 检查文件路径是否有效
+        if not os.path.exists(shp_file_path):
+            feedback.reportError(f"指定的文件路径 {shp_file_path} 不存在!")
+            return {}
+
+        # 加载Shapefile为QgsVectorLayer
+        layer = QgsVectorLayer(shp_file_path, "Vector Layer", 'ogr')
+        print(layer)
+        layer.startEditing()
+        # 获取所有要素(假设网格单元是矩形或正方形)
+        features = list(layer.getFeatures())
+        # 遍历所有网格单元,检查相邻的网格
+        for feature in features:
+            neighbors = self.check_neighbors(feature, layer)
+            for direction, neighbor in neighbors.items():
+                if neighbor:
+                    print(f"")
+                else:
+                    # 更新字段值
+                    feature.setAttribute(direction, "自由")
+                    # 更新要素
+                    layer.updateFeature(feature)
+        # 保存编辑
+        layer.commitChanges()
+        return {
+            "状态": "处理成功"
+        }

+ 264 - 0
processing/tools/Raster/MergeMetadata.py

@@ -0,0 +1,264 @@
+# -*- 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')
+
+
+

+ 151 - 78
processing/tools/Raster/Metadata.py

@@ -8,81 +8,154 @@ __revision__ = '1.0'
 
 from osgeo import gdal, osr
 import openpyxl
-
-# 读取 .img 文件
-filename = 'E:\\projects\\遥感技术部需求\\批量生成元数据\\数据\\分幅成果\\DEM成果\\DEM\\NH48G051053DEMK.img'
-dataset = gdal.Open(filename)
-
-if dataset:
-    # 获取栅格数据的基本信息
-    print(f"Driver: {dataset.GetDriver().ShortName}")
-    print(f"Raster Size: {dataset.RasterXSize} x {dataset.RasterYSize}")
-    print(f"Number of Bands: {dataset.RasterCount}")
-
-    # 获取地理信息
-    geo_transform = dataset.GetGeoTransform()
-    if geo_transform:
-        print(f"GeoTransform: {geo_transform}")
-
-    # 获取投影信息
-    projection = dataset.GetProjection()
-    print(f"Projection: {projection}")
-    # 获取栅格数据的投影信息
-
-    # 创建 SpatialReference 对象并加载投影信息
-    spatial_ref = osr.SpatialReference()
-    spatial_ref.ImportFromWkt(projection)  # 或者使用 ImportFromProj4(projection)
-
-    # 检查投影是否为 UTM 类型,UTM 投影的中央子午线位于 PROJ.4 中的 +lon_0 参数
-    if spatial_ref.IsProjected():
-        proj4 = spatial_ref.ExportToProj4()
-        print(f"Projection PROJ.4: {proj4}")
-
-        # 寻找中央子午线 (lon_0)
-        central_meridian = None
-        for param in proj4.split():
-            if param.startswith('+lon_0'):
-                central_meridian = param.split('=')[1]
-                break
-
-        if central_meridian:
-            print(f"Central Meridian: {central_meridian}")
-        else:
-            print("No central meridian found.")
-
-    # 获取栅格数据的各个波段元数据
-    for band_index in range(1, dataset.RasterCount + 1):
-        band = dataset.GetRasterBand(band_index)
-        print(f"\nBand {band_index}:")
-        print(f"Data type: {gdal.GetDataTypeName(band.DataType)}")
-        print(f"Band min: {band.GetMinimum()}, Band max: {band.GetMaximum()}")
-        print(f"Band NoData Value: {band.GetNoDataValue()}")
-        print(f"Band metadata: {band.GetMetadata()}")
-else:
-    print("Failed to load image.")
-
-
-
-
-
-
-#操作excel
-# filename = 'example.xlsx'
-# workbook = openpyxl.load_workbook(filename)
-# sheet = workbook.active
-# # 定义替换数据的字典
-# replacements = {
-#     "{{name}}": "Bob",
-#     "{{age}}": "40",
-#     "{{address}}": "Paris"
-# }
-# # 遍历单元格,替换占位符
-# 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, replacement)
-# # 保存修改后的文件
-# workbook.save('modified_example.xlsx')
+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')

BIN
processing/tools/Raster/NH48G051053DEMK.xlsx


BIN
processing/tools/Raster/NH48G053055DEMK.xlsx


+ 431 - 0
processing/tools/Raster/makeMetadate.py

@@ -0,0 +1,431 @@
+# -*- 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}个元数据表。"
+        }

BIN
processing/tools/Raster/modified_example.xlsx


BIN
processing/tools/Raster/symbology-style.db


+ 146 - 0
processing/tools/Raster/utils.py

@@ -0,0 +1,146 @@
+# -*- coding: utf-8 -*-
+
+__author__ = 'wanger'
+__description__ = '批量生成栅格数据元数据相关工具'
+__date__ = '2024-11-25'
+__copyright__ = '(C) 2024 by siwei'
+__revision__ = '1.0'
+
+import os
+
+from osgeo import osr
+
+from qgis._core import QgsCoordinateReferenceSystem
+from qgis.core import (
+    QgsRasterLayer,
+    QgsVectorLayer,
+    QgsSpatialIndex,
+    QgsFeatureRequest,
+    QgsGeometry
+)
+
+
+# 通过文件路径来获取文件名和文件后缀
+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 getCentralMeridian(spatial_ref):
+    proj4 = spatial_ref.ExportToProj4()
+    print(f"Projection PROJ.4: {proj4}")
+    central_meridian = None
+    for param in proj4.split():
+        if param.startswith('+lon_0'):
+            central_meridian = param.split('=')[1]
+            break
+    if central_meridian:
+        return central_meridian
+    else:
+        print("No central meridian found.")
+        return None
+
+
+# 通过栅格数据与矢量数据叠加分析  获取相交面积最大的一块  返回要素
+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
+
+# 处理数值
+def processingNumericalValues(value):
+    print(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)