# Python处理PlanetScope影像制作长时序NDVI完整指南 ## 1. PlanetScope数据特性与NDVI计算原理 ### 1.1 PlanetScope影像特点 PlanetScope卫星星座提供高时空分辨率的遥感数据，具有以下关键特征： | 特性 | 参数值 | 说明 | |------|--------|------| | 空间分辨率 | 3-5米 | 高于Landsat和Sentinel-2 | | 重访周期 | 每天 | 可实现密集时间序列分析 | | 波段配置 | 蓝、绿、红、近红外 | 标准四波段多光谱 | | 数据格式 | GeoTIFF | 标准地理空间数据格式 | ### 1.2 NDVI计算公式 归一化差异植被指数（NDVI）是植被监测的核心指标，计算公式为： \[ \text{NDVI} = \frac{\text{NIR} - \text{Red}}{\text{NIR} + \text{Red}} \] 其中NIR代表近红外波段，Red代表红色波段。 ## 2. 数据获取与预处理 ### 2.1 数据下载与组织 ```python import os import glob import numpy as np import rasterio from rasterio.plot import show import matplotlib.pyplot as plt from datetime import datetime # 创建数据目录结构 def create_directory_structure(base_path): """创建标准化的数据目录结构""" directories = { 'raw_data': '原始影像数据', 'processed': '预处理后数据', 'ndvi': 'NDVI计算结果', 'timeseries': '时间序列数据', 'output': '最终输出' } for key, dir_name in directories.items(): path = os.path.join(base_path, dir_name) os.makedirs(path, exist_ok=True) print(f"创建目录: {path}") return directories # 示例使用 base_path = "/path/to/your/planetscope_data" dirs = create_directory_structure(base_path) ``` ### 2.2 数据质量检查与筛选 ```python import json from osgeo import gdal def check_planet_metadata(image_path): """检查PlanetScope影像元数据质量""" try: with rasterio.open(image_path) as src: # 获取基本元数据 metadata = { 'acquisition_date': src.tags().get('ACQUISITION_DATE', 'Unknown'), 'cloud_cover': float(src.tags().get('CLOUD_COVER', 1.0)), 'sun_azimuth': float(src.tags().get('SUN_AZIMUTH', 0)), 'sun_elevation': float(src.tags().get('SUN_ELEVATION', 0)), 'bands': src.count, 'crs': src.crs, 'bounds': src.bounds } # 检查数据有效性 data = src.read(1) valid_pixels = np.count_nonzero(~np.isnan(data)) total_pixels = data.size metadata['valid_ratio'] = valid_pixels / total_pixels return metadata except Exception as e: print(f"读取影像 {image_path} 失败: {e}") return None def filter_images_by_quality(image_dir, max_cloud_cover=0.3, min_valid_ratio=0.8): """基于云覆盖率和有效像素比例筛选影像""" qualified_images = [] for image_file in glob.glob(os.path.join(image_dir, "*.tif")): metadata = check_planet_metadata(image_file) if metadata and (metadata['cloud_cover'] <= max_cloud_cover and metadata['valid_ratio'] >= min_valid_ratio): qualified_images.append({ 'file_path': image_file, 'metadata': metadata }) print(f"合格影像: {os.path.basename(image_file)}, " f"云覆盖率: {metadata['cloud_cover']:.2%}, " f"有效像素: {metadata['valid_ratio']:.2%}") return sorted(qualified_images, key=lambda x: x['metadata']['acquisition_date']) ``` ## 3. NDVI计算核心算法 ### 3.1 基础NDVI计算函数 ```python def calculate_ndvi(red_band, nir_band, output_path=None): """ 计算NDVI指数 参数: red_band: 红色波段数据 (numpy数组或文件路径) nir_band: 近红外波段数据 (numpy数组或文件路径) output_path: 输出文件路径 (可选) 返回: ndvi_array: NDVI计算结果 """ # 处理输入为文件路径的情况 if isinstance(red_band, str): with rasterio.open(red_band) as src: red_data = src.read(1).astype(np.float32) else: red_data = red_band.astype(np.float32) if isinstance(nir_band, str): with rasterio.open(nir_band) as src: nir_data = src.read(1).astype(np.float32) else: nir_data = nir_band.astype(np.float32) # 避免除零错误 denominator = nir_data + red_data denominator[denominator == 0] = np.nan # 计算NDVI ndvi = (nir_data - red_data) / denominator # NDVI值域约束 ndvi = np.clip(ndvi, -1, 1) # 保存结果 if output_path: save_ndvi_raster(ndvi, output_path, red_band if isinstance(red_band, str) else None) return ndvi def save_ndvi_raster(ndvi_data, output_path, reference_image=None): """保存NDVI结果为GeoTIFF文件""" if reference_image: with rasterio.open(reference_image) as src: profile = src.profile profile.update({ 'dtype': rasterio.float32, 'count': 1, 'compress': 'lzw' }) else: profile = { 'driver': 'GTiff', 'dtype': rasterio.float32, 'count': 1, 'crs': 'EPSG:4326', # 默认坐标系 'compress': 'lzw' } with rasterio.open(output_path, 'w', **profile) as dst: dst.write(ndvi_data.astype(np.float32), 1) print(f"NDVI结果已保存: {output_path}") ``` ### 3.2 批量处理与质量控制 ```python def batch_process_ndvi(image_list, output_dir): """批量处理多期影像计算NDVI""" ndvi_results = [] for i, image_info in enumerate(image_list): image_path = image_info['file_path'] acquisition_date = image_info['metadata']['acquisition_date'] print(f"处理第 {i+1}/{len(image_list)} 期影像: {os.path.basename(image_path)}") try: with rasterio.open(image_path) as src: # 假设波段顺序: [蓝, 绿, 红, 近红外] red_band = src.read(3) # 红色波段 nir_band = src.read(4) # 近红外波段 # 生成输出文件名 output_filename = f"NDVI_{acquisition_date}.tif" output_path = os.path.join(output_dir, output_filename) # 计算NDVI ndvi = calculate_ndvi(red_band, nir_band, output_path) ndvi_results.append({ 'date': acquisition_date, 'ndvi_path': output_path, 'ndvi_data': ndvi, 'mean_ndvi': np.nanmean(ndvi), 'max_ndvi': np.nanmax(ndvi), 'min_ndvi': np.nanmin(ndvi) }) except Exception as e: print(f"处理影像 {image_path} 时出错: {e}") continue return ndvi_results ``` ## 4. 时间序列分析与可视化 ### 4.1 时间序列数据构建 ```python import pandas as pd from scipy import stats def build_ndvi_timeseries(ndvi_results, roi_mask=None): """构建NDVI时间序列数据集""" dates = [] mean_ndvi = [] max_ndvi = [] min_ndvi = [] std_ndvi = [] for result in ndvi_results: ndvi_data = result['ndvi_data'] # 如果提供了感兴趣区域掩膜，则只计算该区域 if roi_mask is not None: ndvi_data = ndvi_data[roi_mask] # 计算统计指标 valid_ndvi = ndvi_data[~np.isnan(ndvi_data)] if len(valid_ndvi) > 0: dates.append(result['date']) mean_ndvi.append(np.mean(valid_ndvi)) max_ndvi.append(np.max(valid_ndvi)) min_ndvi.append(np.min(valid_ndvi)) std_ndvi.append(np.std(valid_ndvi)) # 创建时间序列DataFrame ts_df = pd.DataFrame({ 'date': pd.to_datetime(dates), 'mean_ndvi': mean_ndvi, 'max_ndvi': max_ndvi, 'min_ndvi': min_ndvi, 'std_ndvi': std_ndvi }).sort_values('date').reset_index(drop=True) return ts_df def detect_vegetation_phenology(ts_df, window_size=5): """检测植被物候特征""" # 移动平均平滑 ts_df['smooth_ndvi'] = ts_df['mean_ndvi'].rolling( window=window_size, center=True, min_periods=1 ).mean() # 计算物候指标 ts_df['ndvi_trend'] = ts_df['smooth_ndvi'].diff() ts_df['growth_period'] = ts_df['ndvi_trend'] > 0 # 识别关键物候期 max_ndvi_idx = ts_df['smooth_ndvi'].idxmax() min_ndvi_idx = ts_df['smooth_ndvi'].idxmin() phenology_info = { 'peak_ndvi_date': ts_df.loc[max_ndvi_idx, 'date'], 'peak_ndvi_value': ts_df.loc[max_ndvi_idx, 'smooth_ndvi'], 'min_ndvi_date': ts_df.loc[min_ndvi_idx, 'date'], 'min_ndvi_value': ts_df.loc[min_ndvi_idx, 'smooth_ndvi'], 'annual_amplitude': ts_df.loc[max_ndvi_idx, 'smooth_ndvi'] - ts_df.loc[min_ndvi_idx, 'smooth_ndvi'] } return ts_df, phenology_info ``` ### 4.2 高级可视化分析 ```python def plot_ndvi_timeseries(ts_df, phenology_info, output_path=None): """绘制NDVI时间序列图表""" plt.figure(figsize=(12, 8)) # 绘制原始NDVI和平滑NDVI plt.plot(ts_df['date'], ts_df['mean_ndvi'], 'o-', alpha=0.5, label='原始NDVI', color='lightblue') plt.plot(ts_df['date'], ts_df['smooth_ndvi'], '-', linewidth=2, label='平滑NDVI', color='darkblue') # 标记物候关键期 plt.axvline(phenology_info['peak_ndvi_date'], color='red', linestyle='--', alpha=0.7, label=f"NDVI峰值: {phenology_info['peak_ndvi_value']:.3f}") plt.axvline(phenology_info['min_ndvi_date'], color='green', linestyle='--', alpha=0.7, label=f"NDVI谷值: {phenology_info['min_ndvi_value']:.3f}") # 填充生长季 growth_mask = ts_df['growth_period'] == True plt.fill_between(ts_df['date'], ts_df['smooth_ndvi'].min(), ts_df['smooth_ndvi'], where=growth_mask, alpha=0.3, color='green', label='植被生长期') plt.title('PlanetScope NDVI 时间序列分析', fontsize=14, fontweight='bold') plt.xlabel('日期', fontsize=12) plt.ylabel('NDVI', fontsize=12) plt.legend() plt.grid(True, alpha=0.3) plt.xticks(rotation=45) plt.tight_layout() if output_path: plt.savefig(output_path, dpi=300, bbox_inches='tight') print(f"时间序列图已保存: {output_path}") plt.show() def create_ndvi_animation(ndvi_results, output_gif): """创建NDVI变化动画""" import matplotlib.animation as animation fig, ax = plt.subplots(figsize=(10, 8)) def update(frame): ax.clear() ndvi_data = ndvi_results[frame]['ndvi_data'] date = ndvi_results[frame]['date'] im = ax.imshow(ndvi_data, cmap='YlGn', vmin=-0.2, vmax=0.8) ax.set_title(f'NDVI - {date}', fontsize=12) ax.axis('off') # 添加颜色条 if frame == 0: plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04) return [im] anim = animation.FuncAnimation( fig, update, frames=len(ndvi_results), interval=500, blit=False ) anim.save(output_gif, writer='pillow', fps=2) print(f"NDVI动画已保存: {output_gif}") ``` ## 5. 完整工作流示例 ```python def complete_ndvi_workflow(data_directory, output_directory, max_cloud_cover=0.3, start_date=None, end_date=None): """完整的PlanetScope NDVI时间序列分析工作流""" print("=== PlanetScope NDVI时间序列分析工作流开始 ===") # 1. 创建目录结构 dirs = create_directory_structure(output_directory) # 2. 数据质量筛选 print("\n1. 数据质量筛选...") raw_images = glob.glob(os.path.join(data_directory, "*.tif")) qualified_images = filter_images_by_quality( data_directory, max_cloud_cover=max_cloud_cover ) print(f"原始影像数量: {len(raw_images)}") print(f"合格影像数量: {len(qualified_images)}") # 3. 批量计算NDVI print("\n2. 批量计算NDVI...") ndvi_results = batch_process_ndvi( qualified_images, dirs['ndvi'] ) # 4. 构建时间序列 print("\n3. 构建时间序列...") ts_df = build_ndvi_timeseries(ndvi_results) # 5. 物候分析 print("\n4. 物候特征分析...") ts_df, phenology_info = detect_vegetation_phenology(ts_df) # 6. 可视化结果 print("\n5. 生成可视化结果...") timeseries_plot = os.path.join(dirs['output'], 'ndvi_timeseries.png') plot_ndvi_timeseries(ts_df, phenology_info, timeseries_plot) # 7. 生成动画 animation_path = os.path.join(dirs['output'], 'ndvi_animation.gif') create_ndvi_animation(ndvi_results[:10], animation_path) # 限制前10期避免文件过大 # 8. 保存分析报告 report_path = os.path.join(dirs['output'], 'analysis_report.txt') with open(report_path, 'w', encoding='utf-8') as f: f.write("PlanetScope NDVI时间序列分析报告\n") f.write("=" * 50 + "\n\n") f.write(f"分析时间: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}\n") f.write(f"处理影像数量: {len(ndvi_results)}\n") f.write(f"时间范围: {ts_df['date'].min()} 至 {ts_df['date'].max()}\n") f.write(f"NDVI峰值: {phenology_info['peak_ndvi_value']:.3f} " f"(日期: {phenology_info['peak_ndvi_date']})\n") f.write(f"NDVI年振幅: {phenology_info['annual_amplitude']:.3f}\n") f.write(f"平均NDVI: {ts_df['mean_ndvi'].mean():.3f}\n") print("\n=== 工作流完成 ===") print(f"结果保存路径: {output_directory}") return { 'time_series': ts_df, 'phenology': phenology_info, 'ndvi_results': ndvi_results } # 执行完整工作流 if __name__ == "__main__": data_dir = "/path/to/planetscope/images" output_dir = "/path/to/analysis/results" results = complete_ndvi_workflow(data_dir, output_dir) ``` ## 6. 技术要点与最佳实践 ### 6.1 数据处理关键考虑因素 | 技术环节 | 关键考虑 | 解决方案 | |---------|----------|----------| | 数据预处理 | 辐射定标和大气校正 | 使用Planet官方SDK进行辐射归一化 | | 云掩膜 | 云污染影响NDVI精度 | 结合元数据云标记和质量标志波段 | | 时间一致性 | 不同时相影像可比性 | 应用BRDF校正和相对辐射归一化 | | 存储优化 | 大数据量处理效率 | 使用分块处理和内存映射技术 | ### 6.2 性能优化技巧 ```python # 使用分块处理大影像 def process_large_image_chunked(image_path, chunk_size=1024): """分块处理大尺寸PlanetScope影像""" with rasterio.open(image_path) as src: profile = src.profile ndvi_result = np.zeros((src.height, src.width), dtype=np.float32) for ji, window in src.block_windows(1): red_chunk = src.read(3, window=window) nir_chunk = src.read(4, window=window) # 计算当前块的NDVI ndvi_chunk = calculate_ndvi(red_chunk, nir_chunk) ndvi_result[window.row_off:window.row_off+window.height, window.col_off:window.col_off+window.width] = ndvi_chunk return ndvi_result ``` 本方案提供了从数据准备到高级分析的完整技术流程，涵盖了PlanetScope影像处理的核心环节。通过系统化的质量控制、批量处理和可视化分析，用户可以构建高质量的长时序NDVI数据集，为植被动态监测、农业管理和环境变化研究提供可靠的数据支撑[ref_1]。