摘要
准备工作
import rasterioimport numpy as npimport matplotlib.pyplot as pltfrom matplotlib.colors import LinearSegmentedColormapfrom rasterio.warp import transform_boundsfrom pyproj import CRSimport os# 环境配置方法:'''1. 使用conda安装:conda install -c conda-forge rasterioconda install numpy matplotlibconda install -c conda-forge pyproj2. 使用pip安装:pip install rasteriopip install numpypip install matplotlibpip install pyproj'''
处理步骤
def extract_important_info(info):"""提取重要的MTL信息"""important_fields = {'DATE_ACQUIRED': '获取日期','SCENE_CENTER_TIME': '场景中心时间','CLOUD_COVER': '云覆盖率','SUN_ELEVATION': '太阳高度角','SPACECRAFT_ID': '卫星ID','WRS_PATH': 'WRS路径','WRS_ROW': 'WRS行'}result = {}for eng_key, ch_key in important_fields.items():for key in info:if eng_key in key:result[ch_key] = info[key]breakreturn result
def read_image_info(folder_path):"""读取影像的基本信息"""try:# 查找tif文件tif_files = [f for f in os.listdir(folder_path) if f.endswith('.TIF') and not f.endswith('_QA.TIF')]if not tif_files:print("未找到TIF影像文件")return# 读取第一个找到的TIF文件tif_path = os.path.join(folder_path, tif_files[0])with rasterio.open(tif_path) as src:print("\n影像基本信息:")print("文件名:", tif_files[0])print("图像的宽度 (Width):", src.width)print("图像的高度 (Height):", src.height)print("波段数 (Number of bands):", src.count)print("坐标参考系 (CRS):", src.crs)print("分辨率 (Resolution):", src.res)print("边界 (Bounds):", src.bounds)except Exception as e:print(f"读取影像信息时出错: {str(e)}")
def process_folder(input_folder, output_folder, esun_values):"""处理单个文件夹的数据"""# 查找MTL文件mtl_files = [f for f in os.listdir(input_folder) if f.endswith('_MTL.txt')]if not mtl_files:print(f"未在 {input_folder} 中找到MTL文件")return# 读取MTL信息mtl_path = os.path.join(input_folder, mtl_files[0])mtl_info = read_mtl_info(mtl_path)if not mtl_info:return# 计算日地距离date_acquired = datetime.strptime(mtl_info['DATE_ACQUIRED'], '%Y-%m-%d')doy = date_acquired.timetuple().tm_ydayd = 1 - 0.01672 * np.cos(0.9856 * (doy - 4))# 获取太阳高度角sun_elevation = float(mtl_info['SUN_ELEVATION'])# 处理每个波段tif_files = [f for f in os.listdir(input_folder)if f.endswith('.TIF') and 'B' in fand not f.endswith('_QA.TIF')]for tif_file in tif_files:try:band_number = int(tif_file.split('_B')[1].split('.')[0])# 获取定标参数ml_key = f'RADIANCE_MULT_BAND_{band_number}'al_key = f'RADIANCE_ADD_BAND_{band_number}'if ml_key not in mtl_info or al_key not in mtl_info:print(f"波段 {band_number} 缺少定标参数")continueml = float(mtl_info[ml_key])al = float(mtl_info[al_key])# 处理影像input_path = os.path.join(input_folder, tif_file)output_path = os.path.join(output_folder, f'SR_B{band_number}.tif')process_band(input_path, output_path, ml, al, d,esun_values[band_number], sun_elevation)print(f"完成波段 {band_number} 处理: {tif_file}")except Exception as e:print(f"处理波段 {tif_file} 时出错: {str(e)}")
def clip_with_shapefile(input_raster, output_raster, shapefile):"""使用shapefile裁剪栅格数据"""try:# 读取shapefilegdf = gpd.read_file(shapefile)with rasterio.open(input_raster) as src:if gdf.crs != src.crs:gdf = gdf.to_crs(src.crs)# 获取研究区的几何形状shapes = [feature['geometry'] for _, feature in gdf.iterrows()]# 执行裁剪out_image, out_transform = mask(src, shapes, crop=True)out_meta = src.meta.copy()out_meta.update({"driver": "GTiff","height": out_image.shape[1],"width": out_image.shape[2],"transform": out_transform})# 保存裁剪结果with rasterio.open(output_raster, "w", **out_meta) as dest:dest.write(out_image)print(f"裁剪完成: {output_raster}")return Trueexcept Exception as e:print(f"裁剪过程出错: {str(e)}")return False
WRS：201024
WRS：201023
for band, files in band_groups.items():print(f"\n处理波段: {band}")print(f"发现 {len(files)} 个文件")try:src_files = []for path in files:src = rasterio.open(path)src_files.append(src)# 拼接影像mosaic, out_trans = merge(src_files)out_meta = src_files[0].meta.copy()out_meta.update({"driver": "GTiff","height": mosaic.shape[1],"width": mosaic.shape[2],"transform": out_trans})output_filename = f"mosaic_{satellite}_{band}.TIF"output_path = os.path.join(output_folder, output_filename)
band_files = {}for i in range(2, 8):band_file = os.path.join(satellite_folder, f"mosaic_{os.path.basename(satellite_folder)}_B{i}.TIF")if os.path.exists(band_file):band_files[f"B{i}"] = band_fileif not band_files:print(f"在 {satellite_folder} 中未找到波段文件")return# 统计信息print(f"\n{os.path.basename(satellite_folder)} 波段统计信息:")stats_data = []headers = ["波段", "最小值", "最大值", "平均值"]
for idx, (band_name, file_path) in enumerate(band_files.items()):row = idx // 3col = idx % 3with rasterio.open(file_path) as src:band_data = src.read(1, masked=True)valid_data = band_data.compressed()data_range = np.percentile(valid_data, [2, 98])bins = np.linspace(data_range[0], data_range[1], n_bins)cmap = LinearSegmentedColormap.from_list('custom', colors, N=n_bins)colors_list = [cmap(i) for i in np.linspace(0, 1, n_bins)]n, bins, patches = axs[row, col].hist(valid_data, bins=bins,edgecolor='black',linewidth=1)bin_centers = 0.5 * (bins[:-1] + bins[1:])col_val = np.linspace(0, 1, len(bin_centers))for c, p in zip(col_val, patches):plt.setp(p, 'facecolor', cmap(c))axs[row, col].set_title(f'{band_name} 直方图', fontsize=12, fontweight='bold')axs[row, col].set_xlabel('反射率值', fontsize=10)axs[row, col].set_ylabel('像元数量', fontsize=10)
# 绘制每个波段for idx, (band_name, file_path) in enumerate(band_files.items()):row = idx // 3col = idx % 3with rasterio.open(file_path) as src:band = src.read(1, masked=True)band = np.ma.masked_invalid(band)# 将UTM坐标转换为经纬度bounds = src.boundsleft, bottom, right, top = transform_bounds(src.crs, dst_crs,bounds.left, bounds.bottom,bounds.right, bounds.top)# 标准差拉伸stretched_band = standard_deviation_stretch(band)# 绘制图像im = axs[row, col].imshow(stretched_band, cmap=cmap,extent=[left, right, bottom, top])axs[row, col].set_title(f'{band_name}', fontsize=12, fontweight='bold', pad=10)
# 1. 真彩色合成 (4,3,2)rgb_natural = create_composite_image(band_files['B4'],band_files['B3'],band_files['B2'])axs[0, 0].imshow(rgb_natural, extent=[left, right, bottom, top])axs[0, 0].set_title('真彩色合成 (R:B4 G:B3 B:B2)', fontsize=12, fontweight='bold', pad=10)# 2. 标准假彩色 (5,4,3)rgb_false = create_composite_image(band_files['B5'],band_files['B4'],band_files['B3'])axs[0, 1].imshow(rgb_false, extent=[left, right, bottom, top])axs[0, 1].set_title('标准假彩色 (R:B5 G:B4 B:B3)', fontsize=12, fontweight='bold', pad=10)# 3. 城市假彩色 (7,6,4)rgb_urban = create_composite_image(band_files['B7'],band_files['B6'],band_files['B4'])axs[1, 0].imshow(rgb_urban, extent=[left, right, bottom, top])axs[1, 0].set_title('城市假彩色 (R:B7 G:B6 B:B4)', fontsize=12, fontweight='bold', pad=10)# 4. 农业假彩色 (6,5,4)rgb_agriculture = create_composite_image(band_files['B6'],band_files['B5'],band_files['B4'])axs[1, 1].imshow(rgb_agriculture, extent=[left, right, bottom, top])axs[1, 1].set_title('农业假彩色 (R:B6 G:B5 B:B4)', fontsize=12, fontweight='bold', pad=10)
# NDVI (归一化植被指数)ndvi = (nir - red) / (nir + red)# NDBI (归一化建筑指数)ndbi = (swir1 - nir) / (swir1 + nir)# EVI (增强型植被指数)G = 2.5 # 增益系数C1 = 6 # 气溶胶阻抗系数1C2 = 7.5 # 气溶胶阻抗系数2L_evi = 1 # EVI的土壤调节因子evi = G * (nir - red) / (nir + C1 * red - C2 * blue + L_evi)# SAVI (土壤调节植被指数)savi = ((nir - red) * (1 + L)) / (nir + red + L)
结果分析
参考教程
1.Landsat 8 Data Users Handbook | U.S. Geological Survey
2.Rasterio: access to geospatial raster data
3.NumPy 教程| 菜鸟教程
PyGeoSense
一起用遥感探索蓝星风景
