GEE python：使用 Google Earth Engine (GEE) 和 Xarray 库进行基于NOAA数据的归一化植被指数（NDVI）的分析

本文链接：https://blog.youkuaiyun.com/qq_31988139/article/details/146935531

使用 Google Earth Engine 和 Xarray 进行 NDVI 分析的博客

在这篇博客中，我们将探讨如何使用 Google Earth Engine (GEE) 和 Xarray 库进行归一化植被指数（NDVI）的分析。我们将详细介绍数据获取、处理和可视化的步骤。

环境准备

首先，我们需要导入必要的库并进行环境设置：

import ee
import geemap
import xarray as xr

!pip install xee
import xee

ee.Authenticate()
ee.Initialize(project='ee-bqt2000204051')

在这里，我们导入了 ee（Google Earth Engine API）、geemap（用于地图可视化）、xarray（用于处理多维数组数据）和 xee（用于处理 GEE 数据）。接着，我们进行身份验证并初始化 GEE。

创建地图

接下来，我们创建一个地图，并选择卫星影像作为底图：

map = geemap.Map(basemap='SATELLITE')
map

定义感兴趣区域（ROI）

我们可以通过定义一个多边形来指定我们的感兴趣区域（ROI）：

roi = ee.Geometry.Polygon(
        [[[97.41242270469664, 29.124481129731294],
          [97.41242270469664, 20.952364405499313],
          [106.39923911094664, 20.952364405499313],
          [106.39923911094664, 29.124481129731294]]], None)
roi

获取 NDVI 数据

然后，我们从 NOAA 的 VIIRS 数据集中获取 NDVI 数据，并选择 NDVI 和质量保证（QA）波段，过滤时间范围为 2020 年至 2021 年：

ndvi = ee.ImageCollection("NOAA/CDR/VIIRS/NDVI/V1").select('NDVI', 'QA').filterDate('2020', '2021')
ndvi

云掩膜函数

为了去除云和阴影对 NDVI 数据的影响，我们定义了一个云掩膜函数：

def cloud_mask(img):
    index = img.select('NDVI').multiply(0.0001)
    qa = img.select('QA')
    cloud = qa.bitwiseAnd(1 << 1).neq(0)
    shadow = qa.bitwiseAnd(1 << 2).neq(0)
    mask = cloud.Or(shadow).Not()
    return index.updateMask(mask).copyProperties(img, img.propertyNames())

应用云掩膜

我们将云掩膜函数应用于 NDVI 数据集，以获得清晰的 NDVI 数据：

ndvi_masked = ndvi.map(cloud_mask)
ndvi_masked

转换为 Xarray 数据集

接下来，我们将 NDVI 数据转换为 Xarray 数据集，以便进行进一步分析：

ds = xr.open_dataset(ndvi_masked, engine='ee', crs='EPSG:4326', scale=0.05, geometry=roi)
ds

选择特定月份的数据

我们可以选择特定月份的数据，例如 1 月份：

sub = ds.sel(time=ds.time.dt.month == 1)
sub

可视化 NDVI 数据

接下来，我们使用 Xarray 的绘图功能来可视化 1 月份的 NDVI 数据：

sub.NDVI.plot(x='lon', y='lat', col='time', col_wrap=8, robust=True, cmap='RdYlGn')

提取特定点的数据

我们还可以提取特定地点的 NDVI 数据，例如某个经纬度点：

point = ds.sel(lon=49.398480435524995, lat=37.26737308003883, time=ds.time.dt.month == 1, method='nearest')
point.NDVI.plot()

平滑 NDVI 数据

为了平滑 NDVI 数据，我们可以使用滚动窗口平均值：

rolling = ds.rolling(time=10, min_periods=1, center=True).mean()
ds_filled = ds.fillna(rolling)

可视化填充后的 NDVI 数据

接下来，我们将填充后的 NDVI 数据可视化：

sub2 = ds_filled.sel(time=ds_filled.time.dt.month == 1)

import matplotlib.pyplot as plt

sub2.NDVI.plot(x='lon', y='lat', col='time', col_wrap=8, robust=True, cmap='RdYlGn')
plt.savefig('ndvi_filled.png', dpi=360, bbox_inches='tight')

提取填充后的特定点数据

最后，我们提取填充后的特定地点的 NDVI 数据：

point2 = ds_filled.sel(lon=49.398480435524995, lat=37.26737308003883, time=ds.time.dt.month == 1, method='nearest')
point2.NDVI.plot()

总结

通过以上步骤，我们成功地使用 Google Earth Engine 和 Xarray 对 NDVI 数据进行了分析和可视化。这种方法不仅适用于 NDVI 数据，也可以扩展到其他遥感数据的处理和分析中。希望这篇博客能为你提供有用的信息和灵感！如果你有任何问题，请随时联系我。

全部代码

import ee
import geemap
import xarray as xr

!pip install xee
import xee

ee.Authenticate()
ee.Initialize(project = 'ee-bqt2000204051')

map = geemap.Map(basemap = 'SATELLITE')
map
map = geemap.Map(basemap = 'SATELLITE')
map

roi = ee.Geometry.Polygon(
        [[[97.41242270469664, 29.124481129731294],
          [97.41242270469664, 20.952364405499313],
          [106.39923911094664, 20.952364405499313],
          [106.39923911094664, 29.124481129731294]]], None)#;map.draw_last_feature.geometry()
roi

ndvi = ee.ImageCollection("NOAA/CDR/VIIRS/NDVI/V1").select('NDVI', 'QA').filterDate('2020','2021')
ndvi

def cloud_mask(img):
  index = img.select('NDVI').multiply(0.0001)
  qa = img.select('QA')
  cloud = qa.bitwiseAnd(1 << 1).neq(0)
  shadow = qa.bitwiseAnd(1 << 2).neq(0)
  mask = cloud.Or(shadow).Not()
  return index.updateMask(mask).copyProperties(img, img.propertyNames())

ndvi_masked = ndvi.map(cloud_mask)
ndvi_masked

ds = xr.open_dataset(ndvi_masked, engine = 'ee', crs = 'EPSG:4326', scale = 0.05, geometry = roi )

ds

sub = ds.sel(time = ds.time.dt.month == 1)
sub

sub.NDVI.plot(x = 'lon', y = 'lat', col = 'time', col_wrap = 8, robust =True, cmap = 'RdYlGn')

point = ds.sel(lon = 49.398480435524995 , lat = 37.26737308003883, time = ds.time.dt.month == 1, method ='nearest')

point.NDVI.plot()

rolling = ds.rolling(time = 10, min_periods = 1, center = True).mean()

ds_filled = ds.fillna(rolling)

sub2 = ds_filled.sel(time = ds_filled.time.dt.month == 1)

import matplotlib.pyplot as plt

sub2.NDVI.plot(x = 'lon', y = 'lat', col = 'time', col_wrap = 8, robust = True, cmap = 'RdYlGn')
plt.savefig('ndvi_filled.png', dpi = 360, bbox_inches = 'tight')

point2 = ds_filled.sel(lon = 49.398480435524995 , lat = 37.26737308003883, time = ds.time.dt.month == 1, method ='nearest')

point2.NDVI.plot()