geemap学习笔记028：Landsat8计算时间序列NDVI并导出

news2024/11/26 13:45:07

前言

本节则是以Landsat8影像数据为例，进行NDVI时间序列计算，并将得到的时间序列NDVI进行展示并导出。

1 导入库并显示地图

import ee
import geemap
import datetime
import pandas as pd
import os
ee.Initialize()

2 定义时间范围

# 定义日期范围
start_date = '2018-01-01'
end_date = '2019-01-01'
date_range = pd.date_range(start_date, end_date, freq='MS')  #按月生成的时间范围
date_range  #查看日期

3 定义NDVI函数

def NDVI(image): #定义NDVI函数
    nir = image.select("SR_B5")
    red = image.select("SR_B4")
    ndvi = image.expression(
    "float(SR_B5 - SR_B4)/(SR_B5 + SR_B4)",
    {
     "SR_B5": nir,
     "SR_B4": red
    }
    ).rename('NDVI_V2')
    return ndvi

4 加载数据并循环计算NDVI

# 创建一个空列表，用于存储结果图像
result_images = []
# 定义区域
roi = ee.FeatureCollection("projects/xiaoliuk/assets/shape/yantai_qu") #获取assets研究区FeatureCollection
# 循环处理每月数据
for i in range(len(date_range)-1):
    #获取开始日期和结束日期
    start_month = date_range[i].strftime('%Y-%m-%d')
    end_month = date_range[i+1].strftime('%Y-%m-%d')
    
    # 获取影像集合
    dataset = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2') \
        .filterDate(start_month, end_month) \
        .filter(ee.Filter.lt('CLOUD_COVER', 80)) \
        .filterBounds(roi)  #筛选LC08数据时间和经过roi的影像,云量可以自己设定
    
    # 获取中位数影像
    median_image = dataset.median()

    # 计算 NDVI
    # ndvi = median_image.normalizedDifference(['SR_B5', 'SR_B4']) #也可以直接计算NDVI
    ndvi = NDVI(median_image)

    # 将结果图像添加到列表中
    result_images.append(ndvi)

5

visParam = {'min': -0.2,'max': 0.8,
 'palette' : [
    '#d73027',
    '#f46d43',
    '#fdae61',
    '#fee08b',
    '#d9ef8b',
    '#a6d96a',
    '#66bd63',
    '#1a9850',
]}  #可视化参数

# 在地图上显示结果
Map = geemap.Map()
Map.centerObject(roi, zoom=10)
Map.addLayer(result_images[0].clip(roi), visParam, 'NDVI') #选择第一个月份的结果并裁剪roi展示
Map

计算结果

6 与GEE计算的32天的NDVI值进行比较

NDVI_32days = ee.ImageCollection("LANDSAT/LC08/C01/T1_32DAY_NDVI")\
                     .filterDate('2020-01-01', '2020-02-01')\
                     .select('NDVI')   #这是GEE计算的32天的NDVI值，可以进行比较
NDVI_32days

Map.addLayer(NDVI_32days.map(lambda image: image.clip(roi)), visParam, 'NDVI_32')
Map

可以看出两种结果差别并不是很大

7 循环加载单景影像

for i in range(len(result_images)):  #循环加载单景影像
    image_roi = result_images[i].clip(roi) #每月的裁剪影像
    file_name = date_range[i].strftime('%Y-%m')  #没每景影像的名字
    Map.addLayer(image_roi, visParam, file_name) #将其添加到地图中
Map

8 导出影像

out_tif = os.path.join(r"../03_result/") #结果保存路径

for i in range(len(result_images)-2):  #循环加载单景影像
    image_roi = result_images[i].clip(roi)  #每月裁剪的ndvi影像
    file_name = date_range[i].strftime('%Y-%m') 
    
    file_path = os.path.join(out_tif, file_name + '.tif')
    geemap.download_ee_image(
    image=image_roi,
    filename=file_path,
    region=roi.geometry(),
    crs="EPSG:4326",
    scale=30,
    ) #导出影像