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