【RS】GEE(Python)：栅格计算

在遥感影像处理中，栅格计算是一项至关重要的操作。栅格数据代表了地球表面特定范围内的物理量信息，利用栅格计算可以进行多种分析操作，比如计算植被指数、分类、过滤、组合波段，甚至执行复杂的空间分析任务。本篇教程将详细介绍遥感影像栅格计算的多种操作，涵盖多波段组合、指数计算、条件过滤、空间统计、时序分析等进阶功能。

栅格计算概述

栅格数据由行列组成的像元网格表示。每个像元代表特定区域内的某一属性，比如反射率、温度等。栅格计算则是对这些像元执行各种数学运算或逻辑运算，从而提取出有用的信息。

遥感影像通常由多个波段组成，每个波段对应电磁波谱中的不同部分。对不同波段进行组合、计算是栅格计算的常见形式，比如通过不同波段的组合来计算植被指数，或者利用栅格数据来执行分类、聚合等操作。

基础数学运算

在栅格计算中，常见的数学运算包括相加、相减、相乘和相除。这些运算可以用于对影像的不同波段进行组合和计算，从而得出新的分析结果。例如，NDVI、EVI 等常用指数就是通过栅格的波段间的数学运算得出的。

以下是如何使用这些基础数学运算的具体示例：

栅格相加

相加操作可以用于将两个或多个影像的波段加在一起。例如，某些情况下可以将两个波段的像素值相加，用于亮度分析或总值计算。

python">import ee# 创建两个 Landsat 影像波段
band1 = ee.Image('LANDSAT/LC08/C02/T1_L2/LC08_044034_20200114').select('SR_B4')  # 红光波段
band2 = ee.Image('LANDSAT/LC08/C02/T1_L2/LC08_044034_20200114').select('SR_B5')  # 近红外波段# 栅格相加
added_bands = band1.add(band2)# 可视化参数
vis_params = {'min': 0,'max': 50000
}# 显示图像
map = folium.Map(location=[37.9, -122.3], zoom_start=10)
map.add_ee_layer(added_bands, vis_params, 'Added Bands')
map

栅格相减

相减通常用于计算差异，例如，计算某个区域在不同时间点之间的差异（变化检测），或用于计算波段间的差异，例如，NDVI 计算。

python"># 栅格相减
subtracted_bands = band2.subtract(band1)# 可视化参数
vis_params = {'min': 0,'max': 50000
}# 显示图像
map = folium.Map(location=[37.9, -122.3], zoom_start=10)
map.add_ee_layer(subtracted_bands, vis_params, 'Subtracted Bands')
map

栅格相乘

相乘运算通常用于比例调整，或计算波段间的乘积。例如，可以用来放大影像的亮度值，或者计算某些遥感指数时使用。

python"># 栅格相乘
multiplied_bands = band1.multiply(band2)# 可视化参数
vis_params = {'min': 0,'max': 5000000000
}# 显示图像
map = folium.Map(location=[37.9, -122.3], zoom_start=10)
map.add_ee_layer(multiplied_bands, vis_params, 'Multiplied Bands')
map

栅格相除

相除通常用于波段比例计算，例如常见的归一化植被指数（NDVI）的计算公式：

$\frac{(NIR - Red)}{(NIR + Red)}$

python"># 栅格相除
divided_bands = band2.divide(band1)# 可视化参数
vis_params = {'min': 0,'max': 2
}# 显示图像
map = folium.Map(location=[37.9, -122.3], zoom_start=10)
map.add_ee_layer(divided_bands, vis_params, 'Divided Bands')
map

组合运算：NDVI 计算

NDVI 是最常用的遥感指数之一，它结合了近红外波段和红光波段的相减与相加运算。具体公式为：

$\frac{(B5 - B4)}{(B5 + B4)}$

python"># 计算 NDVI
ndvi = band2.subtract(band1).divide(band2.add(band1)).rename('NDVI')# 可视化参数
vis_params = {'min': -1,'max': 1,'palette': ['blue', 'white', 'green']
}# 显示图像
map = folium.Map(location=[37.9, -122.3], zoom_start=10)
map.add_ee_layer(ndvi, vis_params, 'NDVI')
map

条件运算（例如：掩膜处理）

条件运算通常用于创建掩膜，或者对影像进行过滤。例如，使用 lt()、gt() 进行阈值处理。

python"># 掩膜处理，将值大于2000的像素保留，小于2000的像素设为无效值
masked_band = band1.updateMask(band1.gt(2000))# 可视化参数
vis_params = {'min': 0,'max': 3000
}# 显示图像
map = folium.Map(location=[37.9, -122.3], zoom_start=10)
map.add_ee_layer(masked_band, vis_params, 'Masked Band')
map

栅格统计（例如：求平均值）

可以对栅格数据进行统计操作，例如计算某一地区内的平均值、最大值、最小值等。

python"># 定义感兴趣区
region = ee.Geometry.Rectangle([-123.0, 37.0, -122.0, 38.0])# 计算区域内波段的平均值
mean_value = band1.reduceRegion(reducer=ee.Reducer.mean(),geometry=region,scale=30
)print(mean_value.getInfo())

栅格重采样

对于不同分辨率的影像，可以进行重采样操作，以确保统一的像素大小。

python"># 使用最近邻法进行重采样
resampled_image = band1.resample('nearest').reproject(crs='EPSG:4326', scale=500)# 可视化参数
vis_params = {'min': 0,'max': 3000
}# 显示图像
map = folium.Map(location=[37.9, -122.3], zoom_start=10)
map.add_ee_layer(resampled_image, vis_params, 'Resampled Image')
map

栅格计算的常见操作

波段组合与指数计算

在遥感分析中，指数计算是最常见的栅格计算之一。比如，植被指数 (NDVI) 是通过组合近红外和红光波段计算得出的：

$\text{NDVI} = \frac{(NIR - Red)}{(NIR + Red)}$
常见的指数包括：

NDVI (归一化植被指数)：用于衡量植被活力，计算公式为 (NIR - Red) / (NIR + Red)。
EVI (增强植被指数)：一种更精确的植被指数，考虑了大气校正。
NDWI (归一化水体指数)：用于检测水体。
SAVI (土壤调整植被指数)：适用于土壤较为暴露的区域，加入了土壤校正因子。

例：NDVI 计算

python"># 计算 NDVI，B5 代表近红外波段，B4 代表红光波段
ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI')# 可视化 NDVI
ndvi_vis = {'min': -1, 'max': 1, 'palette': ['blue', 'white', 'green']}
map.add_ee_layer(ndvi, ndvi_vis, 'NDVI')

例：EVI 计算
EVI 的计算公式为：

$\text{EVI} = G \times \frac{(NIR - Red)}{(NIR + C1 \times Red - C2 \times Blue + L)}$

python"># 计算 EVI
evi = image.expression('2.5 * ((B5 - B4) / (B5 + 6 * B4 - 7.5 * B2 + 1))',{'B5': image.select('B5'),  # NIR'B4': image.select('B4'),  # Red'B2': image.select('B2')   # Blue}
).rename('EVI')# 可视化 EVI
evi_vis = {'min': -1, 'max': 1, 'palette': ['white', 'blue', 'green']}
map.add_ee_layer(evi, evi_vis, 'EVI')

空间统计与区域聚合

栅格数据的区域聚合是通过计算某个地理区域内像元的统计信息来提取有意义的结果，比如最大值、最小值、平均值等。这些操作通常用于生成空间统计图，分析特定区域内的变化趋势。

python"># 计算指定区域内的 NDVI 平均值
region = ee.Geometry.Polygon([[[-122.5, 37.0], [-122.5, 38.0], [-121.5, 38.0], [-121.5, 37.0]]])ndvi_mean = ndvi.reduceRegion(reducer=ee.Reducer.mean(),geometry=region,scale=30
)print('NDVI Mean:', ndvi_mean.getInfo())

时序分析与栅格叠加

栅格时序分析用于在不同时间点上监测地表特征的变化。可以通过叠加时间序列影像，计算趋势或变化速率。Google Earth Engine 的 ImageCollection 提供了丰富的时序分析工具。

例：计算时间序列的 NDVI 变化

python"># 对影像集合计算 NDVI，并按时间提取 NDVI 序列
def compute_ndvi(image):return image.normalizedDifference(['B5', 'B4']).rename('NDVI')ndvi_collection = landsat_collection.map(compute_ndvi)# 计算时间序列 NDVI 的平均值
ndvi_timeseries = ndvi_collection.reduce(ee.Reducer.mean())

滑动窗口计算与空间滤波

滑动窗口计算是通过移动窗口计算局部统计值，通常用于图像的平滑、降噪、边缘检测等。这种方法对影像中的局部特征进行空间统计分析非常有效。

python"># 使用卷积滤波进行影像平滑
kernel = ee.Kernel.gaussian(radius=3, sigma=1, units='pixels')
smoothed_image = image.convolve(kernel).rename('Smoothed')# 可视化平滑结果
smooth_vis = {'min': 0, 'max': 3000, 'bands': ['B4', 'B3', 'B2']}
map.add_ee_layer(smoothed_image, smooth_vis, 'Smoothed Image')

分类与监督学习

栅格分类是遥感影像分析的关键技术之一，通过将每个像元分配到一个类别，实现对影像中的不同地物类型进行分类。Google Earth Engine 支持多种分类算法，包括监督分类和非监督分类。

例：K-means 非监督分类

python"># 使用 K-means 进行非监督分类
training_data = image.sample(region=region, scale=30, numPixels=5000)
clusterer = ee.Clusterer.wekaKMeans(5).train(training_data)
classified = image.cluster(clusterer)# 可视化分类结果
classification_vis = {'min': 0, 'max': 4, 'palette': ['blue', 'green', 'red', 'yellow', 'purple']}
map.add_ee_layer(classified, classification_vis, 'K-means Classification')

高效处理与导出栅格

计算完成后，通常需要将结果保存或导出，以便进一步分析或用于报告。Google Earth Engine 支持将影像导出为多种格式，包括 GeoTIFF、CSV、Shapefile 等。

python"># 导出 NDVI 结果为 GeoTIFF
export_task = ee.batch.Export.image.toDrive(image=ndvi,description='NDVI_Export',scale=30,region=region
)export_task.start()

进阶技巧与并行计算

对于大规模遥感数据分析，通常需要提高计算效率。Google Earth Engine 提供了并行处理的能力，可以通过 map 函数对整个影像集合并行执行栅格计算。

python"># 并行计算多个时间点的 NDVI
ndvi_series = landsat_collection.map(lambda img: img.normalizedDifference(['B5', 'B4']).rename('NDVI'))

结论

栅格计算是遥感数据分析的核心技术，通过不同的计算方法和空间分析手段，可以对地表特征进行深度挖掘。通过本文的讲解，我们了解了波段组合、指数计算、条件过滤、空间统计和时序分析等多个功能。在实际应用中，可以根据具体的任务需求选择合适的栅格运算方法，从而为科研、监测、决策等提供有力支持。