一、概述：数据预处理部分主要围绕 MODIS ET（蒸散发）数据的格式转换、质量控制、空间裁剪与拼接展开，目的是将原始 HDF 数据处理为研究区（淄博市）的高质量栅格数据。
二、数据预处理整体逻辑：
预处理流程遵循 “原始数据解析→格式转换→质量控制→空间裁剪→区域合并” 的递进式逻辑，逐步将原始遥感数据转化为符合研究需求的标准化数据，具体步骤如下：
1.原始 HDF 数据转 TIFF 格式：将 HDF 文件中的 ET（蒸散发）和 QC（质量控制）子数据集提取并转换为通用的 TIFF 格式，便于后续处理；
2.基于 QC 的质量掩膜：利用 QC 数据过滤 ET 中的无效值（如噪声、低质量数据），保留可靠数据；
3.研究区空间裁剪：根据淄博市边界（GeoJSON 矢量数据），裁剪 TIFF 数据，仅保留研究区内的像元；
4.分区数据合并：将裁剪后的分区数据拼接为完整的淄博市范围数据，形成全市域的标准化栅格数据。
三、各步骤详细逻辑与知识点
数据准备：
# 全局参数设置
base_path = "Zibo_ET_Annual"
boundary_path = os.path.join(base_path, "boundary", "zibo_city.json")
raw_hdf_path = os.path.join(base_path, "raw_hdf")
# 输出结果路径（Excel、PNG）
excel_output = os.path.join(excel_path,"2024_ET_statistics.xlsx")
plot_output_distribution = os.path.join(plots_path,"plots", "2024_ET_annual_distribution.png")
output_tiff_single_path = os.path.join(output_path, "output_tiff_single")# 已裁剪 已掩膜 淄博市各区市
output_tiff_merge_path = os.path.join(output_path, "merged_tiff") # 已掩膜 已裁剪 淄博市合并
geotiff_path = os.path.join(output_path, "et_tiff") # 转换后的ET TIFF
qc_path = os.path.join(output_path, "qc_tiff") # 转换后的QC TIFF
masked_path = os.path.join(output_path, "masked_tiff") # 掩膜后的ET TIFF 未裁剪
# MODIS数据参数
et_band_name = "ET_500m"
qc_band_name = "ET_QC_500m"
invalid_value = 32767 # 无效值标识
qc_valid_bit = 0 # 用于判断有效性的比特位（需参考官方文档）
qc_valid_value = 0 # 该比特位为0时表示有效
scale_factor = 0.1 # ET数据缩放系数（根据MOD16产品文档设置）
# valid_range = (-32767, 32700) # 有效值范围（超出则视为无效）
fill_values = [32767, 32766, 32765, 32764, 32763, 32762, 32761] # 官方定义的无效值
1. HDF 格式转 TIFF 格式（格式标准化）
核心目的：将 MODIS 原始 HDF 数据（分层数据格式）转换为通用的 TIFF 格式，方便后续栅格数据处理（如裁剪、分析）。
具体操作：
（1）用gdal.Open读取 HDF 文件，通过GetSubDatasets获取包含 ET 和 QC 数据的子数据集（HDF 文件通常包含多个子数据层）；
（2）用gdal.Warp将子数据集转换为 TIFF，设置关键参数：
·坐标系：dstSRS="EPSG:4326"（WGS84 地理坐标系，全球通用的经纬度坐标系）；
·数据类型：ET 设为GDT_Float32（浮点型，适合连续值），QC 设为GDT_UInt16（无符号 整数，适合质量标识）；
·效值：srcNodata和dstNodata分别指定原始无效值（32767）和转换后无效值（np.nan）
关键知识点：
（1）HDF 格式特性：Hierarchical Data Format（分层数据格式），常用于存储多波段 / 多参数遥感数据，需通过子数据集路径访问具体数据；
（2）GDAL 库：地理数据抽象库（Geospatial Data Abstraction Library），gdal.Warp是核心函数，支持格式转换、坐标系转换、重投影等；
（3）坐标系：EPSG:4326 对应 WGS84 坐标系，是遥感数据处理中最常用的地理坐标系（经纬度单位）；
（4）数据类型设计：根据数据特性选择（ET 为连续值用浮点型，QC 为标识用整数型），减少存储冗余并保证精度。
代码部分：
# 第一步：将HDF中的ET和QC分别转为单独的TIFF
input_folder_list = os.listdir(raw_hdf_path)
hdf_files = [os.path.join(raw_hdf_path, f) for f in input_folder_list if f.endswith('.hdf')]
if not hdf_files:
print("未找到HDF文件！")
else:
print(f"找到{len(hdf_files)}个HDF文件")
for hdf_file in hdf_files:
try:
hdf_ds = gdal.Open(hdf_file)
if hdf_ds is None:
print(f"无法打开文件: {hdf_file}")
continue
# 获取所有子数据集
sub_datasets = hdf_ds.GetSubDatasets()
# 查找ET和QC波段
et_subdataset = None
qc_subdataset = None
for sub_name, sub_desc in sub_datasets:
if et_band_name in sub_desc:
et_subdataset = sub_name
elif qc_band_name in sub_desc:
qc_subdataset = sub_name
if et_subdataset is None or qc_subdataset is None:
print(f"在文件 {hdf_file} 中未找到所需波段")
continue
# 构建输出TIFF文件路径
base_filename = os.path.splitext(os.path.basename(hdf_file))[0]
et_tiff_path = os.path.join(geotiff_path, f"{base_filename}_et.tif")
qc_tiff_path = os.path.join(qc_path, f"{base_filename}_qc.tif")
# 转换ET波段为TIFF
gdal.Warp(et_tiff_path, et_subdataset,
dstSRS="EPSG:4326",
outputType=gdal.GDT_Float32,
srcNodata=invalid_value,
dstNodata=np.nan)
# 转换QC波段为TIFF
gdal.Warp(qc_tiff_path, qc_subdataset,
dstSRS="EPSG:4326",
outputType=gdal.GDT_UInt16, # QC通常是无符号整数
srcNodata=invalid_value,
dstNodata=invalid_value)
print(f"成功转换: {base_filename} (ET和QC)")
except Exception as e:
print(f"处理文件时出错: {e}")
finally:
if 'hdf_ds' in locals() and hdf_ds is not None:
hdf_ds = None
et_tiff图片：
qc_tiff图片：
2.QC 质量掩膜（数据净化）
核心目的：利用 QC（质量控制）数据过滤 ET 中的低质量 / 无效值，提升数据可靠性。
具体操作：
（1）无效值过滤：先将 ET 数据中官方定义的无效值（fill_values = [32767, 32766,...]）设为np.nan；
（2）QC 比特位校验：通过比特位操作（qc_data & (1 << qc_valid_bit)）判断 QC 数据中标识有效性的比特位（第 0 位），若不符合有效条件（qc_valid_value=0），则将对应 ET 值设为np.nan；
（3）物理范围过滤：将 ET 值超出合理范围（<0或>1500，根据 8 天 ET 的物理意义，1500 对应 150mm/8 天，为极端上限）的像元设为np.nan；
（4）缩放处理：用scale_factor=0.1将原始整数数据转换为实际物理值（原始数据是缩放后的整数，需还原为 mm 单位）。
关键知识点：
（1）QC 数据作用：Quality Control 数据用于标识对应观测值的可靠性（如传感器误差、云污染等），通常通过比特位存储多维度质量信息；
（2）比特位操作：遥感 QC 数据常用比特位编码（1 个整数的不同比特位代表不同质量指标），需用位运算（&、<<）解析（如第 0 位为 0 表示有效）；
（3）无效值处理：用np.nan标记无效值（适合浮点型数据），避免后续计算受异常值干扰；
（4）数据缩放：遥感数据为减少存储量，常将实际值乘以缩放因子后存储为整数（如 ET 实际值 = 原始值 ×0.1）。
# 第二步：使用QC数据掩膜ET数据
et_tiff_files = [os.path.join(geotiff_path, f) for f in os.listdir(geotiff_path) if f.endswith('_et.tif')]
if not et_tiff_files:
print("未找到ET TIFF文件！")
else:
print(f"开始掩膜处理 {len(et_tiff_files)} 个ET TIFF文件")
for et_tiff_file in et_tiff_files:
try:
base_filename = os.path.basename(et_tiff_file).replace('_et.tif', '')
qc_tiff_file = os.path.join(qc_path, f"{base_filename}_qc.tif")
masked_tiff_file = os.path.join(masked_path, f"{base_filename}_masked.tif")
if not os.path.exists(qc_tiff_file):
print(f"找不到对应的QC文件: {qc_tiff_file}")
continue
# 打开ET和QC文件
with rasterio.open(et_tiff_file) as et_src:
et_data = et_src.read(1).astype(np.float32)
profile = et_src.profile
with rasterio.open(qc_tiff_file) as qc_src:
qc_data = qc_src.read(1)
# 1. 过滤官方无效值（必须在缩放前）
for fill_val in fill_values:
et_data[et_data == fill_val] = np.nan
# 2. 应用QC掩膜（不可靠数据设为nan）
bit_mask = (qc_data & (1 << qc_valid_bit)) != (qc_valid_value << qc_valid_bit)
et_data[bit_mask] = np.nan
# 3. 强制过滤物理异常值（核心修正：阈值从3000→1500）
et_data[et_data > 1500] = np.nan # 原始值>1500→对应>150mm/8天（极端值上限）
et_data[et_data < 0] = np.nan # 排除负值
# 4. 应用缩放因子（转换为实际mm）
masked_et_data = et_data * scale_factor
# 5. 验证处理效果
valid_pixels = masked_et_data[~np.isnan(masked_et_data)]
if len(valid_pixels) > 0:
print(f"处理后8天ET范围: {np.min(valid_pixels):.2f} - {np.max(valid_pixels):.2f} mm") # 需在0-150mm
print(f"处理后8天ET均值: {np.mean(valid_pixels):.2f} mm") # 需在30-80mm
else:
print("警告: 处理后无有效像素！")
# 保存掩膜后的文件
with rasterio.open(masked_tiff_file, 'w', **profile) as dst:
dst.write(masked_et_data, 1)
print(f"成功掩膜处理: {masked_tiff_file}")
except Exception as e:
print(f"掩膜处理时出错: {e}")
未裁剪已掩膜的瓦片截图（示例）：
3. 研究区空间裁剪（区域聚焦）
核心目的：根据淄博市边界矢量数据，从大范围 TIFF 中提取研究区（淄博市）内的像元，排除无关区域。
具体操作：
（1）读取 GeoJSON 格式的淄博市边界数据（boundary_path），解析为几何对象（geoms）；
（2）用rasterio.mask.mask函数裁剪 TIFF：
①输入参数：原始 TIFF 数据、边界几何对象（geo）、crop=True（裁剪至边界范围）、all_touched=True（边界触碰的像元均保留）；
②更新元数据：裁剪后影像的尺寸（height/width）和空间变换参数（transform）需重新计算并写入新 TIFF。
关键知识点：
（1）矢量 - 栅格裁剪：通过矢量边界（点、线、面）提取栅格数据中对应区域的像元，是地理数据预处理的核心步骤；
（2）GeoJSON 格式：轻量级矢量数据格式，用于存储地理要素（点、线、面）及其属性，广泛用于 WebGIS 和数据交换；
（3）Rasterio 库：专注于栅格数据处理的 Python 库，mask函数高效实现矢量裁剪，profile属性存储栅格元数据（尺寸、坐标系、数据类型等）；
（4）空间变换参数（transform）：描述栅格像元与地理坐标的映射关系（原点、分辨率、旋转等），裁剪后需重新计算以保证坐标准确性。
# 从文件读取淄博市边界的GeoJSON数据
try:
with open(boundary_path, 'r', encoding='utf-8') as f:
geojson_data = f.read()
print(f"成功加载边界数据: {boundary_path}")
except Exception as e:
print(f"无法加载边界数据: {e}")
# 使用备用的简化边界（如果文件加载失败）
geojson_data = """
{
"type": "FeatureCollection",
"features": [
{
"type": "Feature",
"geometry": {
"type": "Polygon",
"coordinates": [
[
[117.55, 36.83],
[118.23, 36.83],
[118.23, 37.18],
[117.55, 37.18],
[117.55, 36.83]
]
]
},
"properties": {
"name": "淄博市"
}
}
]
}
"""
# 解析GeoJSON数据
geoms = json.loads(geojson_data)
# 获取所有转换后的TIFF文件
tiff_files = [f for f in os.listdir(masked_path) if f.endswith('.tif')]
if not tiff_files:
print("未找到TIFF文件进行裁剪！")
else:
print(f"找到{len(tiff_files)}个TIFF文件，开始裁剪...")
for tiff_file in tiff_files:
input_tiff = os.path.join(masked_path, tiff_file)
output_clipped_tiff = os.path.join(output_tiff_single_path, f"clipped_{tiff_file}")
try:
# 打开栅格文件
with rasterio.open(input_tiff) as src:
# 处理每个GeoJSON要素（如果有多个，分别处理）
for i, feature in enumerate(geoms['features']):
# 提取单个几何对象
geo = [feature['geometry']]
# 执行裁剪操作
out_image, out_transform = mask(
src,
geo,
all_touched=True,
crop=True,
nodata=src.nodata
)
# 获取原始数据的元数据
meta = src.meta
# 更新元数据（尺寸和变换参数）
meta.update({
'driver': 'GTiff',
'height': out_image.shape[1], # 裁剪后的高度
'width': out_image.shape[2], # 裁剪后的宽度
'transform': out_transform,
'nodata': src.nodata
})
# 如果有多个要素，在输出文件名后加序号区分
if len(geoms['features']) > 1:
output_clipped_tiff = os.path.join(
output_tiff_single_path,
f"clipped_{i + 1}_{tiff_file}"
)
# 保存单个要素的裁剪结果
with rasterio.open(output_clipped_tiff, 'w', **meta) as dst:
dst.write(out_image)
print(f"成功裁剪并保存要素 {i + 1}: {output_clipped_tiff}")
except Exception as e:
print(f"裁剪文件 {tiff_file} 时发生错误: {str(e)}")
已裁剪已掩膜分区TIIFF格式图片（示例）：
4.分区数据合并（区域整合）
核心目的：将淄博市各分区的裁剪结果拼接为完整的全市域 TIFF，便于后续分析。
具体操作：
（1）按时间分组：将同一时间的分区裁剪文件（命名含clipped_前缀）归为一组；
（2）用rasterio.merge.merge函数合并：
①输入多个分区 TIFF，自动拼接为覆盖全市的栅格数据；
②更新元数据：合并后影像的尺寸（全市范围）和空间变换参数，确保像元对齐。
关键知识点：
（1）栅格合并：当研究区跨多个原始栅格（如 MODIS 数据分幅）时，需通过合并拼接为完整区域，merge函数自动处理重叠区域的像元对齐；
（2）文件分组逻辑：通过文件名解析（提取原始时间标识）实现同时间数据分组，确保合并的是同一时段的分区数据；
（3）元数据一致性：合并前后需保证坐标系、数据类型一致，否则会导致拼接错位或数据错误。
# 将各个区的 TIFF 文件合并为一个完整的市域 TIFF 文件
# 裁剪后的文件(分区)所在文件夹
clipped_folder = output_tiff_single_path
# 1. 收集所有裁剪后的文件，并按原始文件名分组
file_groups = {}
for filename in os.listdir(clipped_folder):
if filename.startswith("clipped_") and filename.endswith(".tif"):
# 提取原始文件名（去掉"clipped_数字_"前缀）
original_name = "_".join(filename.split("_")[2:])
# 按原始文件名分组
if original_name not in file_groups:
file_groups[original_name] = []
file_groups[original_name].append(os.path.join(clipped_folder, filename))
# 2. 遍历每组文件，合并成全市TIFF
for original_name, file_paths in file_groups.items():
try:
# 打开该组所有分区文件
src_files = [rasterio.open(path) for path in file_paths]
# 合并（拼接成全市范围）
merged_data, merged_transform = merge(src_files)
# 更新元数据（尺寸、坐标系等）
meta = src_files[0].meta.copy()
meta.update({
"height": merged_data.shape[1],
"width": merged_data.shape[2],
"transform": merged_transform
})
# 保存合并后的全市文件到新文件夹（添加"city_"前缀）
merged_filename = f"city_{original_name}"
merged_path = os.path.join(output_tiff_merge_path, merged_filename)
with rasterio.open(merged_path, "w", **meta) as dst:
dst.write(merged_data)
print(f"合并完成：{merged_filename}（保存路径：{merged_path}）")
# 关闭打开的文件
for src in src_files:
src.close()
except Exception as e:
print(f"合并{original_name}时出错：{str(e)}")
已裁剪、已掩膜、淄博市整体TIFF格式图片（示例）：
总结：通过上述步骤，原始 HDF 数据被逐步转化为 “格式标准化、质量可靠、范围聚焦、区域完整” 的淄博市 ET 栅格数据。整个过程体现了遥感数据预处理的典型思路：格式转换→质量控制→空间约束→区域整合。
全部代码汇总：
import os
import re
import json
import numpy as np
import rasterio
from rasterio.merge import merge
from rasterio.mask import mask
from osgeo import gdal
import datetime
import matplotlib.pyplot as plt
import pandas as pd
# 全局参数设置
base_path = "Zibo_ET_Annual"
boundary_path = os.path.join(base_path, "boundary", "zibo_city.json")
raw_hdf_path = os.path.join(base_path, "raw_hdf")
output_path = os.path.join(base_path, "output")
excel_path = os.path.join(output_path, "excel")
plots_path = os.path.join(output_path, "plots")
# 输出结果路径（Excel、PNG）
excel_output = os.path.join(excel_path,"2024_ET_statistics.xlsx")
plot_output_distribution = os.path.join(plots_path,"plots", "2024_ET_annual_distribution.png")
output_tiff_single_path = os.path.join(output_path, "output_tiff_single")# 已裁剪 已掩膜 淄博市各区市
output_tiff_merge_path = os.path.join(output_path, "merged_tiff") # 已掩膜 已裁剪 淄博市合并
geotiff_path = os.path.join(output_path, "et_tiff") # 转换后的ET TIFF
qc_path = os.path.join(output_path, "qc_tiff") # 转换后的QC TIFF
masked_path = os.path.join(output_path, "masked_tiff") # 掩膜后的ET TIFF 未裁剪
# MODIS数据参数
et_band_name = "ET_500m"
qc_band_name = "ET_QC_500m"
invalid_value = 32767 # 无效值标识
qc_valid_bit = 0 # 用于判断有效性的比特位（需参考官方文档）
qc_valid_value = 0 # 该比特位为0时表示有效
scale_factor = 0.1 # ET数据缩放系数（根据MOD16产品文档设置）
# valid_range = (-32767, 32700) # 有效值范围（超出则视为无效）
fill_values = [32767, 32766, 32765, 32764, 32763, 32762, 32761] # 官方定义的无效值
# 将HDF格式转换为TIFF格式
# 创建输出目录
os.makedirs(geotiff_path, exist_ok=True)
os.makedirs(qc_path, exist_ok=True)
os.makedirs(masked_path, exist_ok=True)
os.makedirs(output_tiff_merge_path, exist_ok=True)
os.makedirs(output_tiff_single_path, exist_ok=True)
os.makedirs(excel_path, exist_ok=True)
os.makedirs(plots_path, exist_ok=True)
# 第一步：将HDF中的ET和QC分别转为单独的TIFF
input_folder_list = os.listdir(raw_hdf_path)
hdf_files = [os.path.join(raw_hdf_path, f) for f in input_folder_list if f.endswith('.hdf')]
if not hdf_files:
print("未找到HDF文件！")
else:
print(f"找到{len(hdf_files)}个HDF文件")
for hdf_file in hdf_files:
try:
hdf_ds = gdal.Open(hdf_file)
if hdf_ds is None:
print(f"无法打开文件: {hdf_file}")
continue
# 获取所有子数据集
sub_datasets = hdf_ds.GetSubDatasets()
# 查找ET和QC波段
et_subdataset = None
qc_subdataset = None
for sub_name, sub_desc in sub_datasets:
if et_band_name in sub_desc:
et_subdataset = sub_name
elif qc_band_name in sub_desc:
qc_subdataset = sub_name
if et_subdataset is None or qc_subdataset is None:
print(f"在文件 {hdf_file} 中未找到所需波段")
continue
# 构建输出TIFF文件路径
base_filename = os.path.splitext(os.path.basename(hdf_file))[0]
et_tiff_path = os.path.join(geotiff_path, f"{base_filename}_et.tif")
qc_tiff_path = os.path.join(qc_path, f"{base_filename}_qc.tif")
# 转换ET波段为TIFF
gdal.Warp(et_tiff_path, et_subdataset,
dstSRS="EPSG:4326",
outputType=gdal.GDT_Float32,
srcNodata=invalid_value,
dstNodata=np.nan)
# 转换QC波段为TIFF
gdal.Warp(qc_tiff_path, qc_subdataset,
dstSRS="EPSG:4326",
outputType=gdal.GDT_UInt16, # QC通常是无符号整数
srcNodata=invalid_value,
dstNodata=invalid_value)
print(f"成功转换: {base_filename} (ET和QC)")
except Exception as e:
print(f"处理文件时出错: {e}")
finally:
if 'hdf_ds' in locals() and hdf_ds is not None:
hdf_ds = None
#
# 第二步：使用QC数据掩膜ET数据
et_tiff_files = [os.path.join(geotiff_path, f) for f in os.listdir(geotiff_path) if f.endswith('_et.tif')]
if not et_tiff_files:
print("未找到ET TIFF文件！")
else:
print(f"开始掩膜处理 {len(et_tiff_files)} 个ET TIFF文件")
for et_tiff_file in et_tiff_files:
try:
base_filename = os.path.basename(et_tiff_file).replace('_et.tif', '')
qc_tiff_file = os.path.join(qc_path, f"{base_filename}_qc.tif")
masked_tiff_file = os.path.join(masked_path, f"{base_filename}_masked.tif")
if not os.path.exists(qc_tiff_file):
print(f"找不到对应的QC文件: {qc_tiff_file}")
continue
# 打开ET和QC文件
with rasterio.open(et_tiff_file) as et_src:
et_data = et_src.read(1).astype(np.float32)
profile = et_src.profile
with rasterio.open(qc_tiff_file) as qc_src:
qc_data = qc_src.read(1)
# 1. 过滤官方无效值（必须在缩放前）
for fill_val in fill_values:
et_data[et_data == fill_val] = np.nan
# 2. 应用QC掩膜（不可靠数据设为nan）
bit_mask = (qc_data & (1 << qc_valid_bit)) != (qc_valid_value << qc_valid_bit)
et_data[bit_mask] = np.nan
# 3. 强制过滤物理异常值（核心修正：阈值从3000→1500）
et_data[et_data > 1500] = np.nan # 原始值>1500→对应>150mm/8天（极端值上限）
et_data[et_data < 0] = np.nan # 排除负值
# 4. 应用缩放因子（转换为实际mm）
masked_et_data = et_data * scale_factor
# 5. 验证处理效果
valid_pixels = masked_et_data[~np.isnan(masked_et_data)]
if len(valid_pixels) > 0:
print(f"处理后8天ET范围: {np.min(valid_pixels):.2f} - {np.max(valid_pixels):.2f} mm") # 需在0-150mm
print(f"处理后8天ET均值: {np.mean(valid_pixels):.2f} mm") # 需在30-80mm
else:
print("警告: 处理后无有效像素！")
# 保存掩膜后的文件
with rasterio.open(masked_tiff_file, 'w', **profile) as dst:
dst.write(masked_et_data, 1)
print(f"成功掩膜处理: {masked_tiff_file}")
except Exception as e:
print(f"掩膜处理时出错: {e}")
# 从文件读取淄博市边界的GeoJSON数据
try:
with open(boundary_path, 'r', encoding='utf-8') as f:
geojson_data = f.read()
print(f"成功加载边界数据: {boundary_path}")
except Exception as e:
print(f"无法加载边界数据: {e}")
# 使用备用的简化边界（如果文件加载失败）
geojson_data = """
{
"type": "FeatureCollection",
"features": [
{
"type": "Feature",
"geometry": {
"type": "Polygon",
"coordinates": [
[
[117.55, 36.83],
[118.23, 36.83],
[118.23, 37.18],
[117.55, 37.18],
[117.55, 36.83]
]
]
},
"properties": {
"name": "淄博市"
}
}
]
}
"""
# 解析GeoJSON数据
geoms = json.loads(geojson_data)
# 获取所有转换后的TIFF文件
tiff_files = [f for f in os.listdir(masked_path) if f.endswith('.tif')]
if not tiff_files:
print("未找到TIFF文件进行裁剪！")
else:
print(f"找到{len(tiff_files)}个TIFF文件，开始裁剪...")
for tiff_file in tiff_files:
input_tiff = os.path.join(masked_path, tiff_file)
output_clipped_tiff = os.path.join(output_tiff_single_path, f"clipped_{tiff_file}")
try:
# 打开栅格文件
with rasterio.open(input_tiff) as src:
# 处理每个GeoJSON要素（如果有多个，分别处理）
for i, feature in enumerate(geoms['features']):
# 提取单个几何对象
geo = [feature['geometry']]
# 执行裁剪操作
out_image, out_transform = mask(
src,
geo,
all_touched=True,
crop=True,
nodata=src.nodata
)
# 获取原始数据的元数据
meta = src.meta
# 更新元数据（尺寸和变换参数）
meta.update({
'driver': 'GTiff',
'height': out_image.shape[1], # 裁剪后的高度
'width': out_image.shape[2], # 裁剪后的宽度
'transform': out_transform,
'nodata': src.nodata
})
# 如果有多个要素，在输出文件名后加序号区分
if len(geoms['features']) > 1:
output_clipped_tiff = os.path.join(
output_tiff_single_path,
f"clipped_{i + 1}_{tiff_file}"
)
# 保存单个要素的裁剪结果
with rasterio.open(output_clipped_tiff, 'w', **meta) as dst:
dst.write(out_image)
print(f"成功裁剪并保存要素 {i + 1}: {output_clipped_tiff}")
except Exception as e:
print(f"裁剪文件 {tiff_file} 时发生错误: {str(e)}")
# 将各个区的 TIFF 文件合并为一个完整的市域 TIFF 文件
# 裁剪后的文件(分区)所在文件夹
clipped_folder = output_tiff_single_path
# 1. 收集所有裁剪后的文件，并按原始文件名分组
file_groups = {}
for filename in os.listdir(clipped_folder):
if filename.startswith("clipped_") and filename.endswith(".tif"):
# 提取原始文件名（去掉"clipped_数字_"前缀）
original_name = "_".join(filename.split("_")[2:])
# 按原始文件名分组
if original_name not in file_groups:
file_groups[original_name] = []
file_groups[original_name].append(os.path.join(clipped_folder, filename))
# 2. 遍历每组文件，合并成全市TIFF
for original_name, file_paths in file_groups.items():
try:
# 打开该组所有分区文件
src_files = [rasterio.open(path) for path in file_paths]
# 合并（拼接成全市范围）
merged_data, merged_transform = merge(src_files)
# 更新元数据（尺寸、坐标系等）
meta = src_files[0].meta.copy()
meta.update({
"height": merged_data.shape[1],
"width": merged_data.shape[2],
"transform": merged_transform
})
# 保存合并后的全市文件到新文件夹（添加"city_"前缀）
merged_filename = f"city_{original_name}"
merged_path = os.path.join(output_tiff_merge_path, merged_filename)
with rasterio.open(merged_path, "w", **meta) as dst:
dst.write(merged_data)
print(f"合并完成：{merged_filename}（保存路径：{merged_path}）")
# 关闭打开的文件
for src in src_files:
src.close()
except Exception as e:
print(f"合并{original_name}时出错：{str(e)}")
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
