前言
本节则是以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,
) #导出影像
后记
大家如果有地信遥感方面的问题需要请教或者有项目需要合作,可以在闲 鱼软件的用户中搜索:遥感GIS工作室,请认准头像,谢谢。