如何在R中打开AVIRIS高光谱影像并绘制数据散点图

AVIRIS影像读取可行方法 #

AVIRIS高光谱产品主流存储为两类：带.hdr头文件的ENVI标准二进制格式、NetCDF4格式的新一代AVIRIS-NG产品，对应读取方法如下：

• 通用ENVI格式读取：优先用terra包，内存管理效率高，适配TB级大体积高光谱影像，不会轻易出现内存溢出问题。
  操作代码如下：

  # 首次使用先安装包
  install.packages("terra")
  # 加载包
  library(terra)
  # 读取时直接选中同目录下的.hdr头文件即可，无需手动选中二进制影像文件
  aviris_data <- rast("你的文件存储路径/aviris影像名.hdr")
  # 读取后执行以下命令验证数据，可查看波段数、空间范围、波长元数据
  print(aviris_data)
  

• NetCDF格式AVIRIS-NG产品读取：用ncdf4包直接解析nc文件变量。
  操作代码如下：

  install.packages("ncdf4")
  library(ncdf4)
  # 建立nc文件连接
  nc_con <- nc_open("你的文件存储路径/aviris_ng影像名.nc")
  # 查看文件内存储的变量名，反射率/辐射率数据一般对应"reflectance"/"radiance"字段
  names(nc_con$var)
  # 读取目标变量数据
  ref_data <- ncvar_get(nc_con, "reflectance")
  # 读取完成后关闭连接
  nc_close(nc_con)
  # 后续可根据需要将矩阵转成terra的SpatRaster对象做后续处理
  

读取注意事项：必须保证.hdr头文件和对应的二进制影像文件在同一文件夹下，且文件名前缀完全一致，否则会出现波段识别错误、数值乱码问题。AVIRIS原始产品默认字节顺序为大端（头文件中byte order对应值为1），如果读出来数值全为异常值，优先检查头文件中数据类型、字节顺序参数是否匹配。如果头文件丢失，可直接参考对应AVIRIS产品的官方参数手动编写hdr文本文件即可正常读取。

高光谱数据散点图绘制实现 #

高光谱散点图最常见的场景是两个波段的像元值关联散点、地物光谱特征散点，实现流程如下：

1. 数据预处理
   读取后的栅格数据需要先提取目标波段值，过滤背景空值，为了避免全量像元绘图卡顿，建议先做随机抽样：

   # 获取所有波段对应的中心波长
   wave <- as.numeric(names(aviris_data))
   # 定位需要绘图的两个波段，比如680nm红光波段、800nm近红外波段
   b_red <- which.min(abs(wave - 680))
   b_nir <- which.min(abs(wave - 800))
   # 提取两个波段的像元值，转成绘图用数据框
   plot_df <- data.frame(
     red_value = values(aviris_data[[b_red]], mat = FALSE),
     nir_value = values(aviris_data[[b_nir]], mat = FALSE)
   )
   # 过滤空值、0值背景
   plot_df <- na.omit(plot_df)
   plot_df <- plot_df[plot_df$red_value > 0 & plot_df$nir_value > 0, ]
   # 随机抽取2%的像元用于绘图，大幅提升渲染速度，不影响整体分布趋势
   plot_df <- plot_df[sample(nrow(plot_df), nrow(plot_df)*0.02), ]
   

2. 基础绘图（Base R方案，无需额外装包）

   plot(plot_df$red_value, plot_df$nir_value,
        pch = 16, cex = 0.3, col = rgb(0,0,0,0.15), # 半透明设置解决点重叠遮挡问题
        xlab = "680nm红光波段反射率",
        ylab = "800nm近红外波段反射率",
        main = "AVIRIS影像红光-近红外波段像元散点图")
   # 添加线性拟合线
   abline(lm(nir_value ~ red_value, data = plot_df), col = "firebrick", lwd = 2)
   

3. 更美观的ggplot2方案

   install.packages("ggplot2")
   library(ggplot2)
   ggplot(plot_df, aes(x = red_value, y = nir_value)) +
     geom_point(alpha = 0.1, size = 0.5) +
     geom_smooth(method = "lm", color = "firebrick", linewidth = 1) +
     labs(x = "680nm红光波段反射率",
          y = "800nm近红外波段反射率",
          title = "AVIRIS高光谱影像波段散点图") +
     theme_bw()
   

• 小提示：如果需要绘制单类地物的全波段光谱散点，可以先通过掩膜提取目标地物的所有像元光谱，转成长格式数据框后，将波长映射为x轴、反射率映射为y轴即可。

内容的提问来源于stack exchange，提问作者Aly Hass