从手机拍照到卫星成像：用Python+GDAL带你入门遥感数据处理（附实战代码）
·
从手机拍照到卫星成像：用Python+GDAL带你入门遥感数据处理（附实战代码）
当我们用手机拍摄一张照片时，其实已经在进行最简单的遥感实践——记录物体反射的可见光信息。这种日常行为与专业卫星遥感的核心原理惊人地相似：都是通过传感器捕捉电磁波信息来认识世界。本文将用Python和GDAL库搭建一座桥梁，让开发者能够快速上手处理真实的卫星影像数据。
1. 环境配置与数据准备
1.1 安装GDAL及其Python绑定
GDAL(Geospatial Data Abstraction Library)是处理地理空间数据的瑞士军刀。在Python环境中，我们推荐使用conda进行安装，它能自动解决复杂的依赖关系：
conda create -n rs python=3.8
conda activate rs
conda install -c conda-forge gdal
验证安装是否成功：
from osgeo import gdal
print(gdal.__version__)
提示：Windows用户若遇到路径问题，可将GDAL的bin目录添加到系统环境变量PATH中。
1.2 获取卫星影像数据
Landsat和Sentinel系列卫星提供了免费的开源数据。以Landsat 8为例，可以从USGS EarthExplorer下载：
- 注册USGS EarthExplorer账号
- 选择数据集"Landsat Collection 2 Level-2"
- 划定感兴趣区域和时间范围
- 下载包含所有波段的.tar.gz压缩包
典型Landsat 8文件结构：
LC08_L2SP_118038_20201020_20201105_02_T1/
├── LC08_L2SP_118038_20201020_20201105_02_T1_ANG.txt
├── LC08_L2SP_118038_20201020_20201105_02_T1_MTL.txt
├── LC08_L2SP_118038_20201020_20201105_02_T1_QA_RADSAT.tif
└── Band files (B1-B11)
2. 影像读取与基本操作
2.1 加载多波段影像
使用GDAL读取多波段TIFF文件并转换为numpy数组：
import numpy as np
from osgeo import gdal
def load_bands(band_paths):
bands = []
for path in band_paths:
ds = gdal.Open(path)
bands.append(ds.GetRasterBand(1).ReadAsArray())
ds = None # 显式释放资源
return np.stack(bands, axis=0)
# 示例：加载可见光波段(B2蓝, B3绿, B4红)
band_paths = ['B2.TIF', 'B3.TIF', 'B4.TIF']
bands = load_bands(band_paths)
print(f"波段数组形状：{bands.shape}") # 输出：(3, 行数, 列数)
2.2 生成真彩色合成图像
将红、绿、蓝波段组合显示为自然色彩：
import matplotlib.pyplot as plt
def stretch_8bit(arr, percentile=2):
"""将数据拉伸到0-255范围"""
min_val = np.percentile(arr, percentile)
max_val = np.percentile(arr, 100-percentile)
return np.uint8(255 * (arr - min_val) / (max_val - min_val))
rgb = np.zeros((bands.shape[1], bands.shape[2], 3), dtype=np.uint8)
rgb[..., 0] = stretch_8bit(bands[2]) # 红波段
rgb[..., 1] = stretch_8bit(bands[1]) # 绿波段
rgb[..., 2] = stretch_8bit(bands[0]) # 蓝波段
plt.figure(figsize=(10,10))
plt.imshow(rgb)
plt.axis('off')
plt.title('真彩色合成图像')
plt.show()
注意：卫星数据通常存储为16位整数，直接显示会呈现全黑，必须进行适当拉伸。
3. 植被指数计算与分析
3.1 NDVI计算原理
归一化差异植被指数(NDVI)是最常用的植被健康度指标，计算公式为：
NDVI = (NIR - Red) / (NIR + Red)
其中：
- NIR：近红外波段(Landsat 8的B5)
- Red：红色波段(Landsat 8的B4)
NDVI值域为[-1,1]，典型值范围：
- 水体：<0.1
- 裸土：0.1-0.2
- 植被：0.2-0.8（值越高表示植被越茂盛）
3.2 Python实现NDVI计算
def calculate_ndvi(red_band, nir_band):
"""计算NDVI并处理无效值"""
red = red_band.astype(np.float32)
nir = nir_band.astype(np.float32)
# 避免除以零
denominator = (nir + red)
denominator[denominator == 0] = 1e-10
ndvi = (nir - red) / denominator
# 将无效值设为-1
ndvi[(nir == 0) & (red == 0)] = -1
return ndvi
red = gdal.Open('B4.TIF').GetRasterBand(1).ReadAsArray()
nir = gdal.Open('B5.TIF').GetRasterBand(1).ReadAsArray()
ndvi = calculate_ndvi(red, nir)
# 可视化
plt.figure(figsize=(10,10))
plt.imshow(ndvi, cmap='RdYlGn', vmin=-0.5, vmax=1)
plt.colorbar(label='NDVI值')
plt.title('NDVI分布图')
plt.axis('off')
plt.show()
3.3 进阶应用：植被覆盖度估算
基于NDVI可进一步估算植被覆盖度(FVC)：
def fvc_from_ndvi(ndvi, soil_ndvi=0.2, veg_ndvi=0.86):
"""估算植被覆盖度"""
fvc = (ndvi - soil_ndvi) / (veg_ndvi - soil_ndvi)
fvc = np.clip(fvc, 0, 1) # 限制在0-1范围
return fvc
fvc = fvc_from_ndvi(ndvi)
4. 影像裁剪与保存结果
4.1 基于地理坐标的影像裁剪
使用GDAL的Warp功能实现精确地理裁剪：
from osgeo import gdal, osr
def clip_by_coordinates(input_path, output_path, bounds, epsg=4326):
"""按地理坐标范围裁剪影像"""
# bounds格式：(min_lon, min_lat, max_lon, max_lat)
options = gdal.WarpOptions(
outputBounds=bounds,
outputBoundsSRS=f'EPSG:{epsg}',
resampleAlg='bilinear'
)
gdal.Warp(output_path, input_path, options=options)
# 示例：裁剪杭州市区范围
hangzhou_bounds = (120.0, 30.1, 120.3, 30.4)
clip_by_coordinates('B4.TIF', 'B4_clip.tif', hangzhou_bounds)
4.2 保存计算结果为GeoTIFF
将NDVI结果保存为具有地理参考的TIFF文件：
def save_geotiff(output_path, array, template_path):
"""使用模板文件的投影和转换信息保存数组"""
template_ds = gdal.Open(template_path)
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create(
output_path,
xsize=array.shape[1],
ysize=array.shape[0],
bands=1,
eType=gdal.GDT_Float32
)
out_ds.SetGeoTransform(template_ds.GetGeoTransform())
out_ds.SetProjection(template_ds.GetProjection())
out_ds.GetRasterBand(1).WriteArray(array)
out_ds.FlushCache()
out_ds = None
save_geotiff('NDVI.tif', ndvi, 'B4.TIF')
5. 进阶处理技巧
5.1 批量处理多时相数据
分析同一区域不同时间的植被变化：
import glob
import pandas as pd
def batch_ndvi_analysis(data_folder):
"""批量处理多时相NDVI"""
dates = []
mean_ndvis = []
for scene in glob.glob(f'{data_folder}/LC08_*'):
date = scene.split('_')[3][:8] # 从文件名提取日期
red_path = f'{scene}/B4.TIF'
nir_path = f'{scene}/B5.TIF'
red = gdal.Open(red_path).GetRasterBand(1).ReadAsArray()
nir = gdal.Open(nir_path).GetRasterBand(1).ReadAsArray()
ndvi = calculate_ndvi(red, nir)
# 排除无效值(-1)后计算均值
valid_ndvi = ndvi[ndvi != -1]
if len(valid_ndvi) > 0:
dates.append(pd.to_datetime(date, format='%Y%m%d'))
mean_ndvis.append(np.mean(valid_ndvi))
return pd.DataFrame({'Date': dates, 'NDVI': mean_ndvis}).sort_values('Date')
ndvi_series = batch_ndvi_analysis('Landsat_data')
ndvi_series.plot(x='Date', y='NDVI', title='NDVI时间序列变化')
5.2 使用Rasterio简化操作
Rasterio是基于GDAL的更友好Python接口：
import rasterio
from rasterio.plot import show
# 读取波段
with rasterio.open('B4.TIF') as src:
red = src.read(1)
profile = src.profile # 保存元数据
# 修改元数据后保存NDVI
profile.update(dtype=rasterio.float32, count=1)
with rasterio.open('NDVI_rasterio.tif', 'w', **profile) as dst:
dst.write(ndvi.astype(rasterio.float32), 1)
# 交互式显示
fig, ax = plt.subplots(figsize=(10,10))
show(ndvi, transform=src.transform, ax=ax, cmap='RdYlGn')
plt.colorbar(ax=ax.get_images()[0], label='NDVI值')
6. 常见问题解决方案
6.1 坐标系转换问题
当需要将结果与其他GIS数据叠加时，可能需要进行坐标系转换：
def reproject_raster(input_path, output_path, target_epsg):
"""重投影栅格数据"""
options = gdal.WarpOptions(
dstSRS=f'EPSG:{target_epsg}',
resampleAlg='bilinear'
)
gdal.Warp(output_path, input_path, options=options)
# 示例：转为Web墨卡托投影(EPSG:3857)
reproject_raster('NDVI.tif', 'NDVI_3857.tif', 3857)
6.2 处理大内存影像
对于超大影像，可分块处理以避免内存溢出：
def block_process(input_path, output_path, process_func, block_size=1024):
"""分块处理大影像"""
src_ds = gdal.Open(input_path)
cols = src_ds.RasterXSize
rows = src_ds.RasterYSize
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create(
output_path,
cols,
rows,
1,
gdal.GDT_Float32
)
out_ds.SetGeoTransform(src_ds.GetGeoTransform())
out_ds.SetProjection(src_ds.GetProjection())
for i in range(0, rows, block_size):
for j in range(0, cols, block_size):
# 计算当前块的实际大小
actual_block_x = min(block_size, cols - j)
actual_block_y = min(block_size, rows - i)
# 读取块数据
block = src_ds.GetRasterBand(1).ReadAsArray(
j, i, actual_block_x, actual_block_y
)
# 处理并写入结果
processed_block = process_func(block)
out_ds.GetRasterBand(1).WriteArray(
processed_block, j, i
)
out_ds.FlushCache()
src_ds = None
out_ds = None
# 使用示例：分块计算NDVI
def ndvi_process(block):
"""假设block包含红和近红外波段交替的行"""
red = block[::2]
nir = block[1::2]
return calculate_ndvi(red, nir)
block_process('combined_bands.tif', 'ndvi_blocks.tif', ndvi_process)
掌握这些核心技能后，你已经可以处理大多数基础遥感分析任务。在实际项目中，我通常会先对小范围测试区域验证处理流程，确认无误后再扩展到全图处理，这能有效避免长时间运行后才发现参数问题的情况。
更多推荐
7
0- 0
已为社区贡献34条内容
所有评论(0)