GEE Python中按年掩蔽MODIS NDVI影像集并提取50%阈值日期

需求说明 #

• 处理MODIS日尺度NDVI影像，识别每年NDVI在当年最小值与最大值区间内达到50%阈值对应的儒略日
• 生成包含该50%NDVI值及对应儒略日的年度影像集
• 基于年度阈值影像集，按年份批量掩蔽原始日时间序列影像

已实现的基础代码 #

# 加载MODIS NDVI影像集
collection = ee.ImageCollection('MODIS/MCD43A4_006_NDVI') \
                  .filter(ee.Filter.date('2002-01-01', '2019-06-30'))\
                  .select("NDVI")

# 为影像添加儒略日波段
def addDateBand(image):
  doy = image.date().getRelative('day', 'year')
  doyBand = ee.Image.constant(doy).uint16().rename('doy')
  return image.addBands(doyBand)

collection_withdates = collection.map(addDateBand)

# 计算50%阈值的表达式
calc_exp = '(b("NDVI_min") + (((b("NDVI_max") - b("NDVI_min")) * 0.5)))'

# 生成年度50%阈值影像集
def makeYearlyComposite(date):
  date = ee.Date(date)
  return (collection_withdates \
          .filterDate(date, date.advance(1, 'year')) \
          .reduce(ee.Reducer.minMax())\
          .expression(calc_exp).rename('NDVI_50')) \
          .set("system:index", date.format("YYYY")) 

start = ee.Date('2002-01-01')
end = ee.Date('2019-06-30')
n_yrs = end.difference(start, 'year').subtract(1)
years = ee.List.sequence(0, n_yrs).map(lambda n : start.advance(n, 'year'))
                                           
calc_composite = ee.ImageCollection(years.map(makeYearlyComposite))

批量掩蔽解决方案 #

方法一：直接匹配年份执行掩蔽 #

通过为每幅日影像匹配对应年份的阈值影像，直接执行掩蔽操作：

def maskImageWithYearlyThreshold(image):
    # 获取当前影像的年份
    image_year = image.date().format('YYYY')
    # 筛选对应年份的阈值影像
    threshold_img = calc_composite.filter(ee.Filter.eq('system:index', image_year)).first()
    # 执行掩蔽：保留NDVI大于等于50%阈值的区域，可根据需求调整条件（如lte）
    return image.updateMask(image.select('NDVI').gte(threshold_img))

# 对整个日影像集应用批量掩蔽
masked_collection = collection_withdates.map(maskImageWithYearlyThreshold)

方法二：预生成阈值字典优化效率 #

先将年度阈值影像集转换为年份-影像的字典，减少循环中的筛选操作，提升效率：

# 将年度阈值影像集转为年份与影像对应的字典
year_threshold_features = calc_composite.map(lambda img: 
    ee.Feature(None, {'year': img.get('system:index'), 'image': img})
)
year_threshold_dict = year_threshold_features.reduceColumns(
    ee.Reducer.toList('year'), ee.Reducer.toList('image')
).getInfo()
year_threshold_dict = dict(zip(year_threshold_dict['list'][0], year_threshold_dict['list'][1]))

# 定义掩蔽函数
def maskByYearlyThreshold(image):
    year = image.date().format('YYYY').getInfo()
    threshold_img = ee.Image(year_threshold_dict[year])
    # 执行掩蔽并添加阈值波段（可选）
    return image.updateMask(image.select('NDVI').gte(threshold_img)).addBands(threshold_img)

# 生成批量掩蔽后的影像集
masked_collection = collection_withdates.map(maskByYearlyThreshold)

扩展：添加50%阈值对应的儒略日 #

如果需要在年度影像集中同时保留50%阈值对应的儒略日，可修改makeYearlyComposite函数，找到NDVI最接近阈值的影像并提取其儒略日：

def makeYearlyComposite(date):
    date = ee.Date(date)
    yearly_col = collection_withdates.filterDate(date, date.advance(1, 'year'))
    # 计算年度NDVI的min和max
    min_max = yearly_col.reduce(ee.Reducer.minMax())
    ndvi_50 = min_max.expression(calc_exp).rename('NDVI_50')
    
    # 计算每幅影像NDVI与阈值的差值绝对值，找到最接近阈值的影像
    def calcDiff(img):
        diff = img.select('NDVI').subtract(ndvi_50).abs().rename('diff')
        return img.addBands(diff)
    closest_img = yearly_col.map(calcDiff).sort('diff').first()
    
    # 提取对应儒略日并添加到年度影像
    doy_50 = closest_img.select('doy').rename('doy_50')
    return ndvi_50.addBands(doy_50).set("system:index", date.format("YYYY"))

# 重新生成包含儒略日的年度阈值影像集
calc_composite_with_doy = ee.ImageCollection(years.map(makeYearlyComposite))

内容的提问来源于stack exchange，提问作者M. McPartland