实验要求
使用Landsat合成5-9月生长季均值EVI数据 下载并查看EVI数据结果
任务分解
计算EVI指数，所以需要使用地表反射率数据（surface reflectance，SR）； Landsat数据筛选、去云 计算EVI指数 预览EVI计算结果 下载数据，检查本地的结果
geemap合成EVI数据并下载
启动geemap
import ee
import geemap
import os
ee.Initialize()
设置研究区范围
ROI = geemap.shp_to_ee('I:/geemap/roi/JJJ.shp')\
.geometry()
Landsat去云
现在Landsat数据更新到了Level2 Collection2，相比老版本的数据去云代码略有不同，去云的具体介绍可以参考“GEE学习室”公众号发布的推文：G006-GEE之Landsat影像数据集（最全）去云教程
这里我把JS版GEE的代码转成了geemap版本，供大家调用：
# 去云函数，适用于 Landsat 4、5、7
def rmL457CloudNew(image):
cloudShadowBitMask = (1 << 4)
cloudsBitMask = (1 << 3)
qa = image.select('QA_PIXEL')
mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0).And(qa.bitwiseAnd(cloudsBitMask).eq(0))
return image.updateMask(mask).copyProperties(image, ["system:time_start"])
# 云遮罩函数
def rmL8Cloud(image):
cloudShadowBitMask = (1 << 4)
cloudsBitMask = (1 << 3)
qa = image.select('QA_PIXEL')
mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0).And(qa.bitwiseAnd(cloudsBitMask).eq(0))
return image.updateMask(mask).copyProperties(image, ["system:time_start"])
# 应用缩放因子函数
def applyScaleFactors(image):
opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2)
thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0)
return image.addBands(opticalBands, None, True).addBands(thermalBands, None, True)
Landsat EVI指数计算
增强植被指数（Enhanced Vegetation Index——EVI）简介
计算公式：
特征： EVI常用于LAI值高，即植被茂密区; 值的范围是-1~1，一般绿色植被区的范围是0.2-0.8
增强型植被指数（EVI）算法是遥感专题数据产品中生物物理参数产品中的一个主要算法，可以同时减少来自大气和土壤噪音的影响，稳定地反应了所测地区植被的情况。红光和近红外探测波段的范围设置更窄，不仅提高了对稀疏植被探测的能力，而且减少了水汽的影响，同时，引入了蓝光波段对大气气溶胶的散射和土壤背景进行了矫正。
代码计算
自定义一个EVI计算函数进行计算：
# 计算EVI的函数
def calculateEVI(image):
evi = image.expression(
'2.5 * ((NIR - Red) / (NIR + 6 * Red - 7.5 * Blue + 1))',
{
'NIR': image.select('SR_B5'),
'Red': image.select('SR_B4'),
'Blue': image.select('SR_B2')
}
).rename('EVI')
return image.addBands(evi)
数据合成并预览
选取Landsat8数据 按研究区筛选影像 按时间筛选，这里是2021年5月1日至2021年9月30日，注意，GEE中 filterDate中日期是左闭右开区间去云 应用比例系数 计算EVI 生长季均值合成 按研究区范围裁剪 预览数据，加连续图例
# 处理Landsat-8数据
LT8IC = (
ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
.filterBounds(ROI)
.filterDate('2021-05-01', '2021-10-01')
.map(rmL8Cloud)
.map(applyScaleFactors)
)
LT8_EVI = (
LT8IC
.map(calculateEVI) # 添加EVI波段
.mean()
.clip(ROI)
)
#可视化显示
vis_EVI = {
"min": -1,
"max": 1,
"palette": ["ffffff", "ce7e45", "df923d", "f1b555", "fcd163", "99b718", "74a901",
"66a000", "529400", "3e8601", "207401", "056201", "004c00", "023b01",
"012e01", "011d01", "011301"],
}
#print(LT8_EVI)
# Add layers to the map
Map = geemap.Map(center=(40, 100), zoom=4)
Map.centerObject(ROI)
Map.addLayer(LT8_EVI.select('EVI'), vis_EVI, 'EVI')
Map.addLayer(ROI, {'color': 'red'}, "roi")
# 添加连续图例 (颜色条)
Map.add_colorbar_branca(
colors=vis_EVI["palette"], # 调色板
vmin=vis_EVI["min"], # 最小值
vmax=vis_EVI["max"], # 最大值
label="EVI", # 图例标题
position="bottomright", # 图例位置
)
Map
下载EVI合成结果
geemap可以将数据直接下载到本地，这点我非常喜欢，比网页版的方便。
output_dir = 'I:/geemap/LT8 EVI/'
if not os.path.exists(output_dir):
os.makedirs(output_dir)
# 导出影像到 Google Drive
geemap.download_ee_image(
image=LT8_EVI.select('EVI'),
filename=output_dir+"2021JJJGrowSeason.tif",
region=ROI,
scale=30,
crs="EPSG:4326"
)
前面有Landsat4/5/7的去云函数，把LT8的代码稍微改改就能用于4、5、7的数据处理，这里不再赘述。
本地查看和检查
QGIS查看EVI
把下载好的数据直接拖到QGIS中查看，可以发现，灰白一片，咋回事呢，这是因为有异常值干扰了显示
修改一下符号化效果，把显示的值域修改到（-1，1），改一个看着顺眼的颜色。
R语言查看EVI
使用R语言terra包可以查看EVI下载结果，并绘制直方图
library(terra)
library(tidyverse)
EVI_2021 = rast("./LT8 EVI/2021JJJGrowSeason.tif")
summary(EVI_2021)
hist(EVI_2021)
从直方图可以看出，大部分还是集中在正确的值域范围里的，但是可能存在少量异常值，这就会对QGIS预览结果造成干扰，实际使用时可以自定义显示值域，把异常值的显示干扰去除。
R语言去除异常值
在R语言中，可以写一个自定义函数，将[-1,1]值域以外的数值设为NA，使用栅格计算的方法应用函数。下面是两种去除异常值的函数，供大家参考：
ifelse判断去除异常值
ifelse这种判断语句在R语言中执行的效率很低，执行很慢。虽然能达到效果，但是效率较差。
fun_rm = function(x){
ifelse(x>=-1 & x<=1, x, NA)
}
t1=proc.time()
EVI_2021_rm = app(EVI_2021, fun_rm, cores = 56)
t2=proc.time()
t=t2-t1
print(paste0('执行时间：',t[3][[1]],'秒'))
# writeRaster(EVI_2021_rm, "./EVI.tif")
summary(EVI_2021_rm)
向量化函数去除异常值
在R语言中，使用ifelse函数进行条件判断确实可能会比较慢，尤其是当处理大量数据时。为了提高效率，我们可以使用向量化操作来替代ifelse。向量化操作可以利用R的内部优化，通常比逐个元素的循环或条件判断要快得多。
以下是优化后的代码：
fun_rm2 <- function(x) {
x[x < -1 | x > 1] <- NA
x
}
t1=proc.time()
EVI_2021_rm2 = app(EVI_2021, fun_rm2, cores = 56)
t2=proc.time()
t=t2-t1
print(paste0('执行时间：',t[3][[1]],'秒'))
summary(EVI_2021_rm2)
统计一下计算结果发现，两种方法计算结果一致，都完成了异常值去除工作。
参考文献
G006-GEE之Landsat影像数据集（最全）去云教程 https://blog.sciencenet.cn/blog-3424049-1264006.html
