01 説明

1.1 邏輯和流程

簡要流程：

1. 獲取2024年覆蓋北京奧林匹克森林公園的所有Sentinel-2影像
2. 對所有不同時間段的影像分別計算NDVI
3. 對於同一時間段的影像，取公園內所有像元NDVI值的中位數作為該時間點的NDVI
4. 將所有時間點的NDVI綜合繪製折線圖
5. 地圖上展示公園的真彩色Sentinel-2圖層和NDVI圖層

ps： 提供JS版本和Python版本(在Pycharm中可正常運行, colab未嘗試)參考學習，二者邏輯基本保持一致.

1.2 數據集説明

使用的數據集為: Harmonized Sentinel-2 MSl:MultiSpectral Instrument, Level-2A (SR)(SR表示表面反射率即地表反射率, TOA版本為大氣層頂反射率<包含大氣層的影響>).

使用到的相關波段信息如下：

其中，比例係數0.0001與像元DN值相乘即可得到真正的Sentinel-2表面反射率(原始值通過整數存儲節省存儲,通過0.0001縮放回來)

Sentinel-2中存在屬性CLOUDY_PIXEL_PERCENTAGE表示影像的雲覆蓋率(單位為%)。

02 JS代碼

/*
demo1 利用Sentinel-2監測北京奧林匹克森林公園2024年NDVI變化

北京奧林匹克森林公園經度: 116.388768°E, 緯度: 39.988588°N
*/

// 準備
var start_date = '2024-01-01';
var end_date = '2024-12-31'
var pt = ee.Geometry.Point([116.388768, 39.988588]);  // 定義矢量點
var roi = pt.buffer(2000);  // 2000m緩衝區
var vis_ture_color = {  // 真彩色顯示參數
  bands: ['B4', 'B3', 'B2'],
  min: 0,
  max: 1,
  gamma: 1.6
}
var vis_ndvi = {  // NDVI顯示參數
  bands: 'NDVI',
  min: 0, 
  max: 1,
  palette: ['white', 'green', 'yellow']
}

// 獲取sentinel-2的影像集合
var s2_collection = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
    .filterBounds(roi)  // 按空間範圍篩選
    .filterDate(start_date, end_date)  // 按時間範圍篩選
    .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20));  // 篩選小於20%雲覆蓋的影像
    
// 計算NDVI
var add_ndvi = function(img){
  var ndvi = img.normalizedDifference(['B8', 'B4']).rename('NDVI');  // 計算NDVI
  return img.addBands([ndvi]);  // 返回添加了NDVI波段的img
}
var s2_ndvi = s2_collection.map(add_ndvi).select('NDVI');

// 繪製ndvi
var ndvi_chart = ui.Chart.image.series({
  imageCollection: s2_ndvi,  // 必填
  region: roi,  // 必填
  reducer: ee.Reducer.mean(),  // 默認就是該參數, 可不填
  scale: 30,  // 在30m的分辨率進行reducer(此處即在30m分辨率下進行均值計算), 可不填
  xProperty: 'system:time_start'  // 默認就是該參數
})
// 為繪製的ndvi添加樣式
ndvi_chart.setOptions({
  title: '2024年北京奧林匹克森林公園NDVI時間序列',
  vAxis: {title: 'NDVI(平均值)', viewWindow: {min: 0, max: null}},
  hAxis: {title: '日期', format: 'MM-dd', gridlines: {count: 10}},
  series:{
    0: {
      color: 'green',
      lineWidth: 1,
      pointSize: 1.5,
      pointShape: 'cicrle'
    }
  }
})
// 在控制枱顯示該圖表
print(ndvi_chart)

// 地圖顯示
var img_median = s2_collection.median().clip(roi).multiply(0.0001);  // 中位數合成
var ndvi_median = s2_ndvi.median().clip(roi);  // 中位數合成ndvi
Map.centerObject(roi, 14);
Map.addLayer(img_median, vis_ture_color, '真彩色影像 (median)');
Map.addLayer(ndvi_median, vis_ndvi, 'NDVI (median)');

結果展示：

03 Python代碼

#%% md
# demo1 利用Sentinel-2監測北京奧林匹克森林公園2024年NDVI變化

北京奧林匹克森林公園經度: 116.388768°E, 緯度: 39.988588°N
#%%
import geemap
import ee
from matplotlib import pyplot as plt
from matplotlib import dates as mdates
import pandas as pd
#%%
# 準備
start_date = '2024-01-01'
end_date = '2024-12-31'
pt = ee.Geometry.Point([116.388768, 39.988588])  # 定義矢量點
roi = pt.buffer(2000)  # 做2000m緩衝區視為北京奧林匹克森林公園區域
#%%
# 獲取Sentinel-2影像集合
s2_collection = (ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
                 .filterDate(start_date, end_date)  # 篩選時間範圍
                 .filterBounds(roi)  # 篩選空間範圍
                 .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))  # 保留雲覆蓋率小於20%的影像
                 # .map(lambda img: img.normalizedDifference(['B8', 'B4']).rename('NDVI'))
                 # .select('NDVI')
                 )
#%%
# 計算NDVI
def add_ndvi(img):
    ndvi = img.normalizedDifference(['B8', 'B4']).rename('NDVI')  # 計算NDVI並修改為波段名稱為NDVI
    return img.addBands(ndvi)  # 添加NDVI波段
s2_ndvi = s2_collection.map(add_ndvi).select('NDVI')
#%%
# 獲取roi的NDVI均值(中位數)
def extract_roi_ndvi(img):
    stats = img.reduceRegion(
        reducer=ee.Reducer.median(),  # 計算的統計量
        geometry=roi,  # 統計的區域
        scale=30  # 允許在30m分辨率進行統計而非此處s2的10m分辨率
    )
    
    return ee.Feature(None).set(stats).set('date', img.date().format('YYYY-MM-dd'))

s2_ndvi_f = s2_ndvi.map(extract_roi_ndvi)
ndvi_list = [p['properties'] for p in s2_ndvi_f.getInfo()['features']]
ndvi_df = pd.DataFrame(ndvi_list).sort_values('date').set_index('date')
ndvi_df
#%%
# 繪製NDVI折線圖
plt.rcParams['font.family'] = ['Times New Roman', 'SimSun']  # 可顯示中文(中文宋體-SimSun, 英文新羅馬-Times New Roman)
plt.rcParams['axes.unicode_minus'] = False  # 可顯示負號

fig, ax = plt.subplots(figsize=(12, 6))
ndvi_df.plot(
    ax=ax,
    style='-o',
    color='green',
    title='2024年北京奧林匹克森林公園NDVI時間序列'
)
# 優化XY軸信息
ax.set_xlabel('日期', size=16)
ax.set_ylabel('NDVI(中位數)', size=16)
ax.grid(True)
ax.set_ylim(0, None)
# 優化X軸日期顯示
ax.xaxis.set_major_formatter(mdates.DateFormatter('%m-%d'))
plt.xticks(size=14)

plt.tight_layout()  # 自動調整佈局(緊湊)
plt.show()
#%%
# 地圖Map上顯示NDVI和真彩色底圖
ndvi_median = s2_ndvi.median().clip(roi)
s2_true_color = s2_collection.median().clip(roi).multiply(0.0001)
Map = geemap.Map()
Map.centerObject(roi, 14)  # roi居中顯示
Map.addLayer(ndvi_median, {'min': 0, 'max': 1, 'palette': ['white', 'green', 'yellow']}, 'NDVI (中位數)')
Map.addLayer(s2_true_color, {'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max': 1, 'gamma': 1.6}, 'Sentinel-2 (真彩色)')
#%%
Map

結果展示：

本文由博客一文多發平台 OpenWrite 發佈！