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

# -*- 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)