在metashape中使用反射率标定板为多光谱数据校准

一、光谱数据导入

首先删除多光谱数据中的rgb图像，并把拍摄的多光谱标定板图像也放入文件夹中

打开metashape导入文件夹

添加照片面板中，选择多光谱选项

图像导入完成

二、反射率标定板配置

1.找到反射率面板

打开工具菜单并选择 Calibrate Reflectance 选项（Tools > Calibrate Reflectance...）。按 Locate Panels 按钮，如果未检测出任何标定板信息，则需要手动添加标定板。

2.手动添加遮罩校准图像

工作区中右击images文件夹，点击“添加摄像头组”创建文件夹并将其命名为校准图像

将定标板图像移入新建的校准图像文件夹

使用智能剪刀工具为每张图像创建遮罩，注意创建的遮罩要覆盖除标定板以外的全部区域，因此可以使用智能剪刀框选出标定板所在区域，点击工具栏的“反向选择”（2），再点击“添加选择”（3），即可创建遮罩。

所有标定板图像遮罩添加完毕后，打开工具选择 Calibrate Reflectance 选项（Tools > Calibrate Reflectance...）,点击locate panels,即可看到添加的标定板。

3.输入标定板的反射率

输入标定板各个波段的反射率值（由生产厂家提供），并勾选光传感校准，然后点击ok

4.设置图像亮度

图像亮度设置后，打开 工具 选择 set brightness 选项,点击“预计”然后点击ok即可，至此反射率标定板配置全部完成。

三、数据处理

1.对齐照片

从 “工作流程” 菜单打开 “对齐照片” 对话框（工作流程>对齐照片...），并使用以下参数：

对齐结果如下：

2.优化摄像机

要提高对齐精度，请使用 Tools （工具） 菜单中的 Optimize Cameras  并选择以下参数进行优化 。

3.构建点云

生成点云可以重建更精确的表面，从而提高最终正射镶嵌的质量。

从 工作流程 菜单（工作流>构建点云...）中打开 构建点云 对话框，并使用以下参数：

使用的质量越高，表面就越精确（点数越多），但需要更长的时间来处理。在大多数情况下，中等质量足以进行航空数据处理，尤其是在地形变化较小的情况下。点击ok，处理结果如下：

4.构建 DEM

从 工作流程 菜单（工作流程>构建 DEM...）中打开 Build DEM 对话框，并使用以下参数：

DEM构建结果如下：

5.构建正射镶嵌（Build Orthomosaic）

从 工作流程 菜单（工作流程>>构建正射镶嵌...）中打开 构建正射镶嵌（Build Orthomosaic） 对话框，并使用以下参数：

总尺寸点击 “预计” 按钮即可显示。构建结果如下：

四、计算植被指数

使用 Tools 菜单中的 Set Raster Transform 选项打开 Raster Calculator 对话框（Tools > Set Raster Transform...）

在 Transform 选项卡上，指定要从源数据计算的指数值，本例为计算NDVI。

在Palette面板中点击自动，计算结果如下：

五、处理结果导出

要导出正射镶嵌生成的结果，请使用文件>导出>导出正射镶嵌。

请注意导出正射镶嵌对话框中的栅格变换部分（红框处）以下选项可用：

• 无 - 表示导出的正射镶嵌将包含与源数据对应的相同数量的波段，任何变换公式都将被忽略;
• 索引值（Index value） - 此选项允许保存 栅格计算器 对话框中的变换公式定义的输出波段;
• 索引颜色(Index color) - 根据 Raster Calculator 对话框中的调色板设置将正射镶嵌保存为 RGB 颜色。导出的栅格图像看起来与正射视图模式下的正射镶嵌显示相同，前提是在栅格计算器对话框的变换部分中选中启用选项。

六、使用python查看输出信息

输出时栅格变换(Raster transform)选择了Index value，因此输出为自己设置的指数，即NDVI。使用代码进一步处理tif图像，输出以下信息：
1.图像的波段数、长宽、地理坐标系。
2.每个波段的均值、平均值、中位值、标准差。
3.每个波段在某个值范围内出现最多的值（即计算每个波段值的频率分布）。

"""
本代码可输出以下信息：
1.图像的波段数、长宽、地理坐标系。
2.每个波段的均值、平均值、中位值、标准差。
3.每个波段在某个值范围内出现最多的值（即计算每个波段值的频率分布）。
4.去除图像内的NAN值
"""
import matplotlib.pyplot as plt
from matplotlib import rcParams
# 设置字体为 SimHei (黑体)
rcParams['font.sans-serif'] = ['SimHei']  # 设置中文字体为黑体
rcParams['axes.unicode_minus'] = False  # 解决负号问题

import rasterio
import numpy as np
import matplotlib.pyplot as plt

# 读取多光谱遥感图像
image_path = "C:/Users/Desktop/测试校准/NDVI.tif"  # 替换为图像路径

with rasterio.open(image_path) as src:
    # 获取波段信息
    bands = src.read()  # 读取所有波段数据，返回一个NumPy数组
    band_count = bands.shape[0]  # 波段数目
    width = src.width  # 图像宽度
    height = src.height  # 图像高度
    crs = src.crs  # 获取坐标参考系

    # 获取图像的元数据
    metadata = src.meta
    print("图像元数据:", metadata)
    print(f"波段数: {band_count}")
    print(f"图像长宽: {width}x{height}")
    print(f"地理坐标系: {crs}")

    # 设置绘图的行列数
    rows = int(np.ceil(band_count / 3))  # 每行显示3个波段（可根据需要调整）
    cols = min(band_count, 3)  # 每列最多显示3个波段

    # 创建子图
    plt.figure(figsize=(15, 5 * rows))  # 设置图形大小

    # 计算每个波段的统计信息并绘制图像
    for i in range(band_count):
        band = bands[i]  # 获取第i个波段数据

        # 删除值为65535, -32767和NaN的像素
        band = band[(band != 65535) & (band != -32767) & ~np.isnan(band)]

        # 计算统计值
        mean = np.mean(band)  # 均值
        median = np.median(band)  # 中位值
        std = np.std(band)  # 标准差

        # 计算频率分布
        unique, counts = np.unique(band, return_counts=True)
        most_common_value = unique[np.argmax(counts)]  # 获取出现最多的值

        # 输出每个波段的统计信息
        print(f"\n第{i + 1}波段统计信息:")
        print(f"均值: {mean}")
        print(f"中位值: {median}")
        print(f"标准差: {std}")
        print(f"出现最多的值: {most_common_value}，出现次数: {np.max(counts)}")

        # 自动划分区间：基于分位数
        percentiles = np.percentile(band, [0, 25, 50, 75, 100])  # 计算分位数
        print(f"自动划分的区间（基于分位数）: {percentiles}")

        # 计算各范围内点个数的比例
        # 使用分位数进行区间划分
        bins = percentiles  # 使用分位数作为区间
        counts_in_ranges = np.histogram(band, bins=bins)[0]  # 计算每个区间内的像素个数
        total_count = np.sum(counts_in_ranges)  # 总有效像素个数
        proportions = counts_in_ranges / total_count  # 计算比例

        # 输出各个区间的比例
        print(f"各个范围内点的比例：")
        for j in range(len(bins) - 1):
            print(f"{bins[j]:.2f} - {bins[j + 1]:.2f}: {proportions[j] * 100:.2f}%")

        # 绘制图像
        plt.subplot(rows, cols, i + 1)  # 创建子图
        plt.imshow(bands[i], cmap='gray', vmin=np.min(band), vmax=np.max(band))  # 根据波段的最小值和最大值进行缩放
        plt.colorbar()
        plt.title(f'第{i + 1}波段')

    plt.tight_layout()  # 调整子图间距
    plt.show()

输出如下：