如何用Python读取卫星遥感影像的波段信息并进行植被指数计算？

咱们搞遥感的，平时看图多了都知道，一张普通的RGB影像也就是个“面子工程”，真要想看出点门道，还得看波段。特别是咱们做植被监测、土地利用分类的时候，单看可见光那几个波段是不够的，必须得把近红外（NIR）这些“隐藏属性”挖出来。

今天咱们不整那些虚头巴脑的理论，直接上手用Python搞定最经典的植被指数——NDVI（归一化植被指数）。这可是遥感处理的“Hello World”，学会了这个，什么EVI、SAVI也就是换个公式的事儿。

工欲善其事：数据源与准备

干咱们这行，数据源是个老大难问题。为了演示方便，咱们假设你已经下载好了一份多波段的GeoTIFF影像（比如Landsat或Sentinel-2的数据）。

如果你手头没有现成的数据，可以去USGS或者ESA的官网下载，或者找我要点测试数据也行。这里咱们默认你的影像已经做好了基本的预处理（比如解压、波段叠加），存成了一个多波段的.tif文件。

Python环境里，核心库就用rasterio。这玩意儿比GDAL的Python接口更“Pythonic”，读写操作符合咱们程序员的直觉，处理坐标系统（CRS）也很稳。

核心操作：波段读取

拿到一张影像，第一件事就是看看它里面到底有几个波段，然后把咱们需要的红光波段和近红外波段“抠”出来。

通常情况下，Landsat 8的红光波段是第4波段，近红外是第5波段；而Sentinel-2的红光是B4，近红外是B8。具体是第几层，得看你原始数据怎么存的，千万别读错了，不然算出来的结果就是一堆乱码。

下面这段代码，是咱们读取波段的常规操作：

import rasterio

import numpy as np

# 打开影像文件

tif_path = 'your_image.tif'

with rasterio.open(tif_path) as src:

# 查看波段数，心里有个底

print(f"波段总数: {src.count}")

# 读取红光波段（假设是第3波段）

# 注意：rasterio的索引是从1开始的，不是0

red_band = src.read(3)

# 读取近红外波段（假设是第4波段）

nir_band = src.read(4)

# 获取原始影像的元数据，后面输出结果要用

profile = src.profile

提个醒：读取出来的数据其实是numpy数组。如果你打印出来看全是黑的，别慌，那是数值太小，显示器显示不出来，数据本身是没问题的。

算法实现：NDVI计算

波段读出来了，剩下的就是小学数学题了。NDVI的公式咱们都倒背如流：

公式简单，但编程实现有个坑：分母为0。在遥感影像里，像元值为0的情况很常见（比如背景值）。如果不处理，程序会报错或者出现无效值。

这时候，咱们得利用numpy的掩膜功能，或者简单的浮点数转换来规避。

# 强制转换为浮点型，防止整数除法截断

red = red_band.astype(float)

nir = nir_band.astype(float)

# 允许除法中出现除零警告，但我们要用np.seterr来忽略它，或者后续处理nan

np.seterr(divide='ignore', invalid='ignore')

# 计算NDVI

ndvi = (nir - red) / (nir + red)

# 将无效值（分母为0的地方）设为nan或特定的nodata值

ndvi[np.isnan(ndvi)] = -9999

计算完的ndvi数组，值域理论上在-1到1之间。负值通常表示水体，0左右是裸土，正值越大植被越茂盛。

实例演示：从计算到导出

光算出来看一眼可不行，咱们得把结果存成一个新的GeoTIFF，方便在ArcGIS或QGIS里出图。

这里我把完整的流程串起来，大家可以直接拿去跑：

import rasterio

import numpy as np

# 1. 读取数据

input_path = 'sentinel2_test.tif' # 替换成你的文件名

with rasterio.open(input_path) as src:

# 假设红光是第3波段，近红外是第4波段

red = src.read(3).astype(float)

nir = src.read(4).astype(float)

# 复制元数据

profile = src.profile

# 2. 计算NDVI

np.seterr(divide='ignore', invalid='ignore')

ndvi = (nir - red) / (nir + red)

# 3. 更新元数据并保存

# NDVI是浮点型，我们要更新dtype

profile.update(

dtype=rasterio.float32,

count=1, # 输出只有1个波段

nodata=-9999

)

# 写入新文件

with rasterio.open('ndvi_result.tif', 'w', **profile) as dst:

dst.write(ndvi.astype(rasterio.float32), 1)

print("NDVI计算完成，结果已保存为 ndvi_result.tif")

跑完这段代码，你会在同目录下得到一个ndvi_result.tif。拖进ArcGIS里，套个绿到红的色带，植被覆盖情况一目了然。

碎碎念（注意事项）

虽然代码跑通了，但实际生产中，这几个坑你大概率会踩：

大气校正问题：如果你直接用DN值（原始灰度值）算NDVI，做定性分析（比如看哪块绿、哪块不绿）没问题。但如果要搞定量反演，比如算具体的叶面积指数，必须先做辐射定标和大气校正，把DN值转成反射率。没做校正的DN值算出来的NDVI，数值偏高，千万别拿去发论文。

分辨率统一：Sentinel-2这种数据，红光波段是10米分辨率，但某些植被红边波段是20米。硬算的话，数组形状不匹配程序直接崩。记得重采样到统一分辨率。

内存溢出：要是处理全省甚至全国的镶嵌影像，直接read()会把内存撑爆。这时候得用rasterio的窗口读取，分块计算，这是进阶话题了，以后有机会再聊。

总结

用Python处理遥感影像，其实就是把GIS软件里的“工具箱”变成了几行代码。虽然刚开始学写代码比点鼠标累，但一旦流程通了，成百上千张图的批处理也就是几秒钟的事儿。

把你自己的影像数据带入上述流程，导出 NDVI 即可。除了 NDVI，不妨试试把公式换成 EVI 或 SAVI，看看不同指数对植被信息的提取效果有什么区别。动手试试，比看十遍教程都管用。

更多相关技术内容咨询欢迎前往并持续关注好学星城论坛了解详情。

想高效系统的学习Python编程语言，推荐大家关注一个微信公众号：Python编程学习圈。每天分享行业资讯、技术干货供大家阅读，关注即可免费领取整套Python入门到进阶的学习资料以及教程，感兴趣的小伙伴赶紧行动起来吧。

发表于 2天前
阅读 ( 4 )
分类：Python开发

如何用Python读取卫星遥感影像的波段信息并进行植被指数计算？

你可能感兴趣的文章

相关问题

0 条评论

作家榜 »