遥感技术作为我们提供了从空中俯瞰地球的能力,它通过卫星、飞机等平台获取地表信息,广泛应用于农业监测、环境评估、城市规划等领域。本文将带你了解遥感的基本概念,并通过 Python 代码实践遥感图像处理的核心技术。
什么是遥感?
- 卫星遥感:如 Landsat、Sentinel 系列卫星
- 航空遥感:飞机、无人机搭载传感器
- 地面遥感:高塔、三脚架等固定平台
Python 遥感处理工具
rasterio
:读取和写入遥感图像numpy
:数值计算基础matplotlib
:图像可视化scikit-image
:图像处理算法earthpy
:简化遥感数据分析流程
实践:遥感图像处理基础
import rasterio
import numpy as np
import matplotlib.pyplot as plt
from rasterio.plot import show
import os
# 设置中文显示
plt.rcParams["font.family"] = ["SimHei", "WenQuanYi Micro Hei", "Heiti TC"]
plt.rcParams['axes.unicode_minus'] = False # 解决负号显示问题
def read_remote_sensing_image(file_path):
"""读取遥感图像文件"""
try:
with rasterio.open(file_path) as src:
# 打印图像信息
print(f"图像文件: {file_path}")
print(f"图像宽度: {src.width} 像素")
print(f"图像高度: {src.height} 像素")
print(f"波段数量: {src.count}")
print(f"空间参考: {src.crs}")
print(f"边界范围: {src.bounds}")
# 读取所有波段
img = src.read()
# 转换形状为 (height, width, bands)
img = np.transpose(img, (1, 2, 0))
return img, src
except Exception as e:
print(f"读取图像出错: {e}")
return None, None
def display_rgb_image(img, src, bands=[3,2,1], title="RGB合成图像"):
"""显示RGB合成图像"""
# 确保波段索引从0开始
bands = [b-1 for b in bands]
# 提取RGB波段
rgb_img = img[:, :, bands]
# 归一化到0-1范围以便正确显示
rgb_normalized = (rgb_img - np.min(rgb_img)) / (np.max(rgb_img) - np.min(rgb_img))
plt.figure(figsize=(10, 8))
plt.imshow(rgb_normalized)
plt.title(title)
plt.axis('off')
plt.show()
return rgb_normalized
def calculate_ndvi(img, nir_band=4, red_band=3):
"""计算归一化植被指数NDVI"""
# 确保波段索引从0开始
nir_idx = nir_band - 1
red_idx = red_band - 1
# 提取近红外和红波段
nir = img[:, :, nir_idx].astype(float)
red = img[:, :, red_idx].astype(float)
# 计算NDVI,避免除零错误
ndvi = np.where(
(nir + red) == 0,
0,
(nir - red) / (nir + red)
)
# 显示NDVI图像
plt.figure(figsize=(10, 8))
plt.imshow(ndvi, cmap='Greens')
plt.colorbar(label='NDVI值')
plt.title('归一化植被指数(NDVI)')
plt.axis('off')
plt.show()
return ndvi
def calculate_ndwi(img, green_band=2, nir_band=4):
"""计算归一化水体指数NDWI"""
# 确保波段索引从0开始
green_idx = green_band - 1
nir_idx = nir_band - 1
# 提取绿波段和近红外波段
green = img[:, :, green_idx].astype(float)
nir = img[:, :, nir_idx].astype(float)
# 计算NDWI,避免除零错误
ndwi = np.where(
(green + nir) == 0,
0,
(green - nir) / (green + nir)
)
# 显示NDWI图像
plt.figure(figsize=(10, 8))
plt.imshow(ndwi, cmap='Blues')
plt.colorbar(label='NDWI值')
plt.title('归一化水体指数(NDWI)')
plt.axis('off')
plt.show()
return ndwi
def classify_land_cover(ndvi, ndwi):
"""基于NDVI和NDWI进行简单土地覆盖分类"""
# 创建分类图像
classification = np.zeros_like(ndvi, dtype=np.uint8)
# 水体:NDWI > 0.1
classification[ndwi > 0.1] = 1
# 植被:NDVI > 0.3 且 不是水体
vegetation_mask = (ndvi > 0.3) & (classification != 1)
classification[vegetation_mask] = 2
# 裸地:NDVI <= 0.3 且 不是水体
bare_ground_mask = (ndvi <= 0.3) & (classification != 1)
classification[bare_ground_mask] = 3
# 显示分类结果
plt.figure(figsize=(10, 8))
cmap = plt.cm.get_cmap('Set3', 4)
im = plt.imshow(classification, cmap=cmap)
plt.colorbar(im, ticks=[0, 1, 2, 3])
plt.clim(-0.5, 3.5)
plt.title('土地覆盖分类')
plt.axis('off')
# 添加图例
from matplotlib.patches import Patch
legend_elements = [
Patch(facecolor=cmap(0), label='其他'),
Patch(facecolor=cmap(1), label='水体'),
Patch(facecolor=cmap(2), label='植被'),
Patch(facecolor=cmap(3), label='裸地')
]
plt.legend(handles=legend_elements, loc='lower right')
plt.show()
return classification
def main():
# 这里使用Landsat 8样例数据,实际使用时替换为你的图像路径
# 可以从USGS Earth Explorer下载免费的Landsat数据
image_path = "LC08_L1TP_123032_20200515_20200527_01_T1_B4.TIF"
# 检查文件是否存在,如果不存在则使用多个波段文件
if not os.path.exists(image_path):
print("找不到单个图像文件,尝试使用多波段文件模式")
# 实际应用中,这里应该是你的多波段文件列表
return
# 读取遥感图像
img, src = read_remote_sensing_image(image_path)
if img is None:
return
# 显示真彩色合成图像 (红、绿、蓝波段)
rgb_img = display_rgb_image(img, src, bands=[4,3,2], title="真彩色合成图像")
# 计算并显示NDVI
ndvi = calculate_ndvi(img)
# 计算并显示NDWI
ndwi = calculate_ndwi(img)
# 进行简单的土地覆盖分类
classification = classify_land_cover(ndvi, ndwi)
print("遥感图像处理完成!")
if __name__ == "__main__":
main()
代码解析
获取遥感数据
- USGS Earth Explorer:提供 Landsat 系列、Sentinel-2 等多种卫星数据
- Sentinel Hub:专注于欧空局 Sentinel 系列卫星数据
- NASA Earthdata:包含多种 NASA 卫星数据
- 地理空间数据云:国内的遥感数据共享平台
遥感技术的应用场景
进阶学习方向
- 高级图像分类算法(如随机森林、深度学习)
- 时间序列分析(监测长期变化趋势)
- 三维遥感(激光雷达数据处理)
- 定量遥感(反演地表参数)