使用geetools绘制Landsat NDVI与EVI时序图空白问题排查

我用Google Earth Engine生成Landsat 8的NDVI与EVI时序图表，切换到Landsat 8影像集后图表空白，但MODIS影像集运行正常。已经确认Landsat计算出的NDVI、EVI值存在且合理，求排查问题。

#--------------------------------------------------------------------------------------
#Python Script to Explore Google Earth Engine
#--------------------------------------------------------------------------------------
import ee
import geetools
import pandas as pd
import geopandas as gpd
import sys 
import os 
import matplotlib.pyplot as plt
from datetime import datetime as d
#--------------------------------------------------------------------------------------
#Functions
#--------------------------------------------------------------------------------------
def apply_scale_factors_ls(image):
    #--These are scale factors for landsat
    try:
      optical_bands = image.select('SR_B.*').multiply(0.0000275).add(-0.2)
      image = image.addBands(optical_bands, None, True)
    except:
      pass 
    return image

def apply_scale_factors_modis(image):
    #--These are scale factors for landsat
    try:
      ndvi = image.select('NDVI').multiply(0.0001)
      evi = image.select('EVI').multiply(0.0001)
      image = image.addBands(ndvi, None, True)
      image = image.addBands(evi, None, True)
    except:
      pass 
    return image
  
def calculate_ndvi(img):
    ndvi = img.normalizedDifference(['SR_B5', 'SR_B4']).rename('NDVI')
    #img.addBands(ndvi)
    #img.select().addBands(ndvi)
    return img.addBands(ndvi)

def calculate_evi(img):
    evi_img = img.select(['SR_B5', 'SR_B4', 'SR_B2'], ['nir', 'red', 'blue'])
    evi_img = evi_img.expression(
        '(2.5 * (NIR - RED)) / (NIR + 6 * RED - 7.5 * BLUE + 1)',
        {
            'NIR': evi_img.select('nir'),
            'RED': evi_img.select('red'),
            'BLUE': evi_img.select('blue')
        }
    ).rename('EVI')
    return img.addBands(evi_img)
#----------------------------------------------------------------------------------
#--Current Working code
#----------------------------------------------------------------------------------
ee.Authenticate()
ee.Initialize(project='bold-camera-445217-d4')
#----------------------------------------------
#--Get our imagery collection 
date_start =  "2018-01-01"
date_end =    "2024-12-01"

bc_extent = ee.Geometry.Rectangle([-139.2,47.5,-113.84,60.1]) 
#bc = geojson_eefc(bc_bounds)
#--MODIS Collection 
sensor = 'MODIS/061/MOD13Q1'
ic = (
    #ee.ImageCollection('MODIS/061/MOD13A1') 
    ee.ImageCollection(sensor) 
    .filter(ee.Filter.date(date_start, date_end)) 
    .select(['NDVI', 'EVI'])
    .map(apply_scale_factors_modis)
    .filterBounds(bc_extent)
    )
maxPixels = 10000000
scale = 500    
#---Landsat Collection 
sensor = 'LANDSAT/LC08/C02/T1_L2'
ic = (
    ee.ImageCollection(sensor)
    .filterDate(date_start, date_end)
    .select(['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5']) 
    .map(apply_scale_factors_ls)
    .map(calculate_ndvi)
    .map(calculate_evi)
    .filterMetadata('CLOUD_COVER', 'less_than', 10)
    .filterBounds(bc_extent)
    )
maxPixels = 300000000
scale = 30
#----------------------------------------------
aoi = ee.Geometry.Polygon([[-124.1753755674451, 51.358079942703306], [-124.0019949017465, 51.358079942703306], [-124.0019949017465, 51.45478364277703], [-124.1753755674451, 51.45478364277703], [-124.1753755674451, 51.358079942703306]])

fig, ax = plt.subplots(figsize=(10, 4))
#--Plot the Current AOS Feature   
ic.geetools.plot_dates_by_bands(
  region = aoi,
  reducer = "mean",
  scale = scale, 
  bands = ["NDVI", "EVI"],
  colors = ["#0193c2","#028352"],
  ax = ax,
  dateProperty = "system:time_start",
  maxPixels = maxPixels
)
ax.set_ylabel("Vegetation Index")
ax.set_title(f"Average Vegetation Index Values in AOS Feature \n {sensor}")

plt.show()

对比截图 #

• MODIS NDVI图表正常显示：
• Landsat 8 NDVI图表空白：

问题排查与修复方案 #

核心问题是Landsat影像集的过滤顺序错误，导致无效影像（不在AOI内或云量过高）被提前处理，同时缺少时序排序，影响geetools的图表渲染逻辑。

修复后的Landsat影像集代码段

#---Landsat Collection 
sensor = 'LANDSAT/LC08/C02/T1_L2'
ic = (
    ee.ImageCollection(sensor)
    .filterDate(date_start, date_end)
    .filterBounds(bc_extent)  # 先过滤AOI范围，减少无效计算
    .filterMetadata('CLOUD_COVER', 'less_than', 10)  # 先过滤云量，保留有效影像
    .select(['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5']) 
    .map(apply_scale_factors_ls)
    .map(calculate_ndvi)
    .map(calculate_evi)
    .sort('system:time_start')  # 按时间排序，确保时序逻辑正确
)
maxPixels = 300000000
scale = 30

备选手动绘图方案（规避geetools潜在兼容问题）

如果上述调整后仍有问题，可以手动提取时序数据到Pandas再绘图，更可控：

# 提取AOI内的均值时间序列
def extract_ts(image):
    mean_val = image.reduceRegion(
        reducer=ee.Reducer.mean(),
        geometry=aoi,
        scale=scale,
        maxPixels=maxPixels
    )
    return ee.Feature(None, {
        'NDVI': mean_val.get('NDVI'),
        'EVI': mean_val.get('EVI'),
        'date': ee.Date(image.get('system:time_start')).format('YYYY-MM-dd')
    })

# 获取特征集并转为Pandas DataFrame
ts_features = ic.map(extract_ts).getInfo()['features']
ts_data = pd.DataFrame([feat['properties'] for feat in ts_features])
ts_data['date'] = pd.to_datetime(ts_data['date'])

# 绘制图表
fig, ax = plt.subplots(figsize=(10, 4))
ts_data.plot(x='date', y='NDVI', color='#0193c2', ax=ax, label='NDVI')
ts_data.plot(x='date', y='EVI', color='#028352', ax=ax, label='EVI')
ax.set_ylabel("植被指数")
ax.set_title(f"AOI内平均植被指数时序图 \n {sensor}")
ax.legend()
plt.show()

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