geemap学习笔记027：遥感影像指数计算，并依据研究区进行裁剪

前言

GEE最常用的一个功能就是指数的计算，并且可以进行较大范围的计算，然后根据感兴趣区进行裁剪。常见的指数有归一化差异植被指数 (NDVI)、增强植被指数 (EVI)、叶面积植被指数（LAI）、归一化差值水体指数(NDWI)、归一化建筑指数（NDBI）等等，这篇文章中有各种指数的介绍。本文则是以Landsat 8影像为基础，以NDVI指数为例进行示例计算。

1 导入库并显示地图

import ee
import geemap

ee.Initialize()
Map = geemap.Map()
Map

2 加载Landsat 8影像和感兴趣

# 应用尺度缩放因子
def apply_scale_factors(image):
  optical_bands = image.select('SR_B.').multiply(0.0000275).add(-0.2)
  thermal_bands = image.select('ST_B.*').multiply(0.00341802).add(149.0)
  return image.addBands(optical_bands, None, True).addBands(
      thermal_bands, None, True
  )

roi = ee.FeatureCollection("projects/xiaoliuk/assets/shape/yantai_qu") #获取assets研究区FeatureCollection

dataset = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2') \
    .filterDate('2022-05-01', '2022-10-01') \
    .filter(ee.Filter.lt('CLOUD_COVER', 15)) \
    .filterBounds(roi)  #筛选LC08数据时间和经过roi的影像

dataset = dataset.map(apply_scale_factors) #对数据应用缩放因子

visualization = {
    'bands': ['SR_B4', 'SR_B3', 'SR_B2'],
    'min': 0.0,
    'max': 0.3,
}  

Map = geemap.Map()
Map.centerObject(roi, 8)
Map.add_layer(dataset, visualization, 'original')
Map

3 定义NDVI函数

这里定义了两种NDVI函数，第一种是集成在GEE中，可以直接调用，但并不是每一种指数都集成到了GEE中；第二种是一种自定义的函数，比较直观，可以根据自己的需求进行修改。

def NDVI_V1(image): #定义第一种NDVI函数
    ndvi = image.normalizedDifference(['SR_B5', 'SR_B4']).rename('NDVI_V1')
    return image.addBands(ndvi)

def NDVI_V2(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 image.addBands(ndvi)

4 计算NDVI并且判断两种方法结果是否一致

withNDVI = dataset.map(NDVI_V1).map(NDVI_V2) #对所有的影像进行两种NDVI方法计算

Map.addLayer(withNDVI.select('NDVI_V1','NDVI_V2'), {}, 'ndvi') #将其添加到地图中，对比两种方法的结果是否一致
Map

在图中随便选择了一个点，可以看出，两种计算方式的结果是一致的

5 影像中值合成并且根据感兴趣区进行裁剪

筛选出来的影像是一个ImageCollection，并不是一张图像，因此可以先进行一步合成。

median = dataset.median()  #对图像进行中值合成，成为一幅影像
ndvi = NDVI_V1(median) #利用第一种方法计算NDVI

Map.add_layer(median, visualization, 'median') #显示中值合成的结果
Map

visParam = {'min': -0.2,'max': 0.8,
 'palette' : [
    '#d73027',
    '#f46d43',
    '#fdae61',
    '#fee08b',
    '#d9ef8b',
    '#a6d96a',
    '#66bd63',
    '#1a9850',
]}  #可视化参数
Map.addLayer(ndvi.select('NDVI_V1').clip(roi), visParam, 'ndvi_median') #将其添加到地图中
Map