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