原文：
towardsdatascience.com/downscaling-a-satellite-thermal-image-from-1000-m-to-10-m-python-3b2ed19ff103
作者可视化的从 1000 m 降尺度到 10 m 的 Sentinel-3 热图像。
目录
-
🌅 引言
-
💾 下载 Sentinel-3 (1000 m) 和 Sentinel-2 图像 (10 m)
-
⚙️ Sentinel-3 图像处理
-
🌡️ 温度-NDVI 空间
-
📐 锐化热图像（1000 m 到 10 m）
-
🗺️ 锐化热图像的可视化
-
📄 结论
-
📚 参考文献
🌅 引言
由于提供热图像的卫星在空间和时间分辨率之间的权衡，卫星热图像的降尺度已经被广泛研究。例如，Landsat-8 的重访周期为 16 天，原始热分辨率为 100 米。相比之下，Sentinel-3 可以提供每日热图像，但空间分辨率为 1000 米。
空间和时间分辨率的权衡，图片由作者提供
解决热图像粗分辨率的一种方法可能是发射更多配备热传感器的卫星，例如 2021 年 9 月发射的 NASA 的 Landsat-9。在这种情况下，Landsat-8 和 Landsat-9 的时空分辨率均为 8 天（而不是一颗卫星的 16 天），假设天气晴朗。
然而，正如你可以猜到的，这种方法需要数百万美元的投资和几年的努力。相反，研究人员已经专注于统计方法，将卫星的可见光/近红外（VNIR）波段与更高空间分辨率（但较低时间分辨率）的卫星热图像相关联，以卫星的低空间分辨率（但较高时间分辨率）的热图像。例如，研究表明，从 Sentinel-2 的 VNIR 波段（10m，每 5 天）计算出的归一化植被指数（NDVI）可以与 Sentinel-3 的热图像（1000m，每日）呈负相关。
但我们如何利用这种关系来锐化热图像呢？
简短的回答是，如果 Sentinel-2 的 NDVI 与 Sentinel-3 的热图像之间的相关性足够强，我们就可以将这个方程应用于 10 米分辨率的 NDVI，从而生成 10 米分辨率的热图像。
如果你想了解更多关于卫星波段和传感器光谱的信息，请参考以下内容：
在这篇帖子中，我们将从 Sentinel-3 下载低空间分辨率的红外图像，并从 Sentinel-2 下载高空间分辨率的 VNIR 图像。这两张图像分别由各自的卫星同时捕获。接下来，我们将使用 VNIR 波段（(NIR-红)/(NIR+红)）计算 NDVI，将其提升到 1000 米分辨率，并探索 NDVI 与热波段（两者均为 1000 米分辨率）之间的相关性。最后，我们将使用这种相关性生成 10 米分辨率的温度图。
如果您喜欢开始的部分，继续阅读还有更多内容可以探索。
💾 下载 Sentinel-3（1000 米）和 Sentinel-2 图像（10 米）
我已经写了三篇关于在 R 和 Python 中下载 Sentinel-2 图像，以及使用 Python 下载 Sentinel-3 图像的帖子。我不想在这里重复这些步骤；相反，我建议您参考以下帖子：
使用 R 下载 Sentinel-2 图像：
使用 Python 下载 Sentinel-2 图像：
使用 Python 下载 Sentinel-3 图像：
如果您不想在编写代码时下载带有文字的图片，请查看这篇帖子：
将统计方法应用于降尺度热图像的关键步骤是找到由卫星同时捕获的清晰图像。在下载图像之前，您可以根据日期、云层覆盖和您感兴趣的区域（AOI）过滤元数据。如果您想了解更多关于过滤和与 Sentinel 元数据（AOI、云层覆盖等）一起工作的信息，请阅读以下内容：
在这篇帖子中，我的 AOI 位于加利福尼亚州，我找到了 2023 年 6 月 19 日由 Sentinel-2 和 Sentinel-3 捕获的清晰图像。您可以自由搜索其他位置或日期。但如果您想使用我下载的图像，以下是相关信息：
Sentinel-2: _S2B_MSIL2A_20230620T183919_N0509_R070_T10SFG20230620T224951
卫星 = “SENTINEL-2”
级别 = “S2MSI2A”
研究区域 = “POLYGON ((-121.0616 37.6391, -120.966 37.6391, -120.966 37.6987, -121.0616 37.6987, -121.0616 37.6391))”
_ 开始日期 = “2023–06–19” ; 结束日期 = “2023–06–21”*
Sentinel-3: _S3A_SL_2_LST____20230620T181455_20230620T181755_20230620T201625_0179_100_141_2340_PS1_O_NR004
卫星 = “SENTINEL-3”
级别 = “LST”
研究区域 = “POLYGON ((-121.0616 37.6391, -120.966 37.6391, -120.966 37.6987, -121.0616 37.6987, -121.0616 37.6391))”
_ 开始日期 = “2023–06–19” ; 结束日期 = “2023–06–21”*
在从 Sentinel-2 下载 10 米（我们需要计算 NDVI 的波段）和 Sentinel-3 的热图像后，你应该在你的目录中有这三个文件：
目录快照，图片由作者提供
⚙️ Sentinel-3 图像处理
Sentinel-3 图像覆盖的场景比 Sentinel-2 大得多。由于我们需要 Sentinel-3 每个像素的平均 NDVI 值，我们需要根据 Sentinel-2 图像的范围裁剪 Sentinel-3 图像。为此，第一步是将 Sentinel-3 图像重新投影到与 Sentinel-2 图像相同的投影。这可以通过以下方式完成：
import numpy as np
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling
from pyproj import Transformer
input_raster = 'Sentinel-3_L2_LST_reproj.tif'
output_raster = 'Sentinel-3_L2_LST_reproj_32610.tif'
dst_crs = 'EPSG:32610'
# Read the input raster file
with rasterio.open(input_raster) as src:
# Get the transform, width, and height for the destination CRS
transform, width, height = calculate_default_transform(src.crs, dst_crs, src.width, src.height, *src.bounds)
# Set up the destination
kwargs = src.meta.copy()
kwargs.update({
'crs': dst_crs,
'transform': transform,
'width': width,
'height': height,
'dtype': np.float32,
})
# Create the destination and write the reprojected data
with rasterio.open(output_raster, 'w', **kwargs) as dst:
# Perform the reprojection
for i in range(1, src.count + 1):
reproject(
source=src.read(1).astype(np.float32)* src.scales[0] + src.offsets[0],
destination=rasterio.band(dst, i),
src_transform=src.transform,
src_crs=src.crs,
dst_transform=transform,
dst_crs=dst_crs,
resampling=Resampling.bilinear)
在这个脚本中，我们读取从上一步下载的 Sentinel-3 栅格图像，计算重新投影到指定 CRS 的转换参数，然后导出带有重新投影数据的新的栅格文件。完成这些步骤后，你应该在你的目录中有这四个文件：
目录快照，图片由作者提供
现在 Sentinel-3 热图像已经重新投影到 Sentinel-2 坐标系中，我们可以继续根据 Sentinel-2 数据范围裁剪热图像，以下代码如下：
import rasterio
import numpy as np
# Open the two raster files
with rasterio.open('T10SFG_20230620T183919_B08_10m.jp2') as small_raster:
with rasterio.open('Sentinel-3_L2_LST_reproj_32610.tif') as big_raster:
# Get the extent of the smaller raster
min_x, min_y, max_x, max_y = small_raster.bounds
# Read the data from the bigger raster within the extent of the smaller raster
window = rasterio.windows.from_bounds(min_x, min_y, max_x, max_y, big_raster.transform)
data = big_raster.read(window=window)
# Update the metadata of the bigger raster to match the extent of the smaller raster
clipped_meta = big_raster.meta.copy()
clipped_meta.update({
'height': window.height,
'width': window.width,
'transform': rasterio.windows.transform(window, big_raster.transform),
'dtype': data.dtype
})
# Write the clipped data
with rasterio.open('Sentinel-3_L2_LST_reproj_32610_clipped.tif', 'w', **clipped_meta) as clipped_raster:
clipped_raster.write(data)
在这个脚本中，我们读取两个栅格文件（Sentinle-3 和 Sentinel-2 图像），提取较小栅格的范围（Sentinel-2），从较大栅格（Sentinle-3 热图像）中读取相应的数据，更新裁剪的 Sentinle-3 的元数据以匹配 Sentinel-2 图像，然后将裁剪的热图像写入新的 TIFF 文件。
目录快照，图片由作者提供
让我们并排绘制 NDVI 和裁剪的温度图：
import numpy as np
import matplotlib.pyplot as plt
import rasterio
from rasterio.plot import show
# File paths
red_path = '/content/T10SFG_20230620T183919_B04_10m.jp2'
nir_path = '/content/T10SFG_20230620T183919_B08_10m.jp2'
clipped_temperature_path = '/content/Sentinel-3_L2_LST_reproj_32610_clipped.tif'
# Read raster data
with rasterio.open(red_path) as red_src:
red = red_src.read(1)
with rasterio.open(nir_path) as nir_src:
nir = nir_src.read(1)
with rasterio.open(clipped_temperature_path) as clipped_temp_ds:
clipped_temperature = clipped_temp_ds.read(1)
# Calculate NDVI
ndvi = (nir - red) / (nir + red)
# Plot NDVI and temperature side by side
fig, (ax1,ax2) = plt.subplots(ncols=2, figsize=(12, 6))
# Plot NDVI
im1 = ax1.imshow(ndvi, cmap=ndvi_cmap, vmin=0, vmax=0.6)
ax1.set_title('NDVI', fontweight='bold', fontsize=14)
fig.colorbar(im1, ax=ax1, shrink=0.5)
# Plot clipped temperature
im2= ax2.imshow(clipped_temperature, cmap=ndvi_cmap.reversed(), vmin=300, vmax=315)
ax2.set_title('Clipped Temperature', fontweight='bold', fontsize=14)
fig.colorbar(im2, ax=ax2, shrink=0.5)
plt.show()
Sentinel-2 10 米分辨率的 NDVI（左）和 1000 米分辨率的裁剪 Sentinel-3 热图像（右），由作者可视化
如您所见，我们能够成功根据哨兵-2 地图（本例中的 NDVI）的范围裁剪哨兵-3 热图像。然而，哨兵-2 图像也被截断，如果我们运行区域统计以获得热像素的平均 NDVI 值，我们将得到许多 NaN 值：
10 米（背景）的哨兵-2 NDVI，哨兵-3 像素（网格，1000 米），作者图片
我们将在下一步中通过排除数据框中的 NaN 值来解决这个问题。
🌡️ 温度-NDVI 空间
通过两个清晰的图像，一个是来自哨兵-3 的热图像，另一个是来自哨兵-2 的 VNIR 图像，在同一投影和范围下，我们可以执行区域统计，以获得每个温度像素的 NDVI 平均值。基本上，通过区域统计方法，我们将 10 米的 NDVI 地图聚合到 1000 米，从而在温度-NDVI 空间中绘制 NDVI 值与热值。
如介绍中所述，热值应与 NDVI 值成反比。较高的 NDVI 表示植被覆盖度更大，对应较冷的像素，而较低的 NDVI 表示植被较少或裸土，对应较暖的像素。让我们进行区域统计，以探索我们 AOI 中热值和 NDVI 值之间的空间：
import rasterio
import rasterio.features
import rasterio.mask
import pandas as pd
import geopandas as gpd
import rasterstats
import rasterio
from rasterio.features import shapes
mask = None
# Open the input raster
with rasterio.open('Sentinel-3_L2_LST_reproj_32610_clipped.tif') as src:
# Read the raster band
image = src.read(1).astype(np.float32)* src.scales[0] + src.offsets[0]
results = (
{'properties': {'Temperature': v}, 'geometry': s}
for i, (s, v)
in enumerate(
shapes(image, mask=mask, transform=src.transform)))
geoms = list(results)
gpd_polygonized_raster = gpd.GeoDataFrame.from_features(geoms)
# Open the rasters
with rasterio.open('T10SFG_20230620T183919_B08_10m.jp2') as nir_src:
with rasterio.open('T10SFG_20230620T183919_B04_10m.jp2') as red_src:
# Read the data as float32
nir = nir_src.read(1).astype(np.float32)* nir_src.scales[0] + nir_src.offsets[0]
red = red_src.read(1).astype(np.float32)* red_src.scales[0] + red_src.offsets[0]
# Calculate NDVI
ndvi = (nir - red) / (nir + red)
# Calculate the zonal statistics for each polygon
stats = rasterstats.zonal_stats(gpd_polygonized_raster.geometry, ndvi, affine=nir_src.transform, stats='mean')
# Add the mean value of NDVI to the dataframe
gpd_polygonized_raster['NDVI'] = [s['mean'] for s in stats]
# Save the polygon layer to a shapefile
gpd_polygonized_raster.to_file('output_polygons.shp')
# Create a pandas dataframe from the geodataframe
stats_df = pd.DataFrame(gpd_polygonized_raster.drop(columns='geometry'))
# Print the dataframe
print(stats_df)
在此脚本中，我们将热值栅格多边形化，从哨兵-2 影像中计算 NDVI，并为每个温度像素计算区域统计（平均 NDVI）。
哨兵-3 热值与哨兵-2 NDVI 值对比
如表格所示，并如前所述，由于哨兵-2 图像的截断，一些像素存在 NaN 值，其中我们拥有来自哨兵-3 的热数据。我们将省略包含 NaN 值的哨兵-2 的行：
# Drop rows with NaN values
df_clean = stats_df.dropna(subset=['NDVI'])
df_clean
输出结果将包括：
移除 nan 值后的哨兵-3 热值与哨兵-2 NDVI 值对比
使用这个干净的数据框，我们可以通过以下方式绘制温度-NDVI 空间：
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import linregress
# Create a scatter plot
plt.figure(figsize=(8, 6))
sns.scatterplot(x='NDVI', y='Temperature', data=df_clean, palette='coolwarm')
# Set plot title and axis labels
plt.title('1-1 Plot of NDVI vs Temperature', fontsize=16, fontweight='bold')
plt.xlabel('NDVI', fontsize=14)
plt.ylabel('Temperature', fontsize=14)
# Show the plot
plt.show()
输出结果将包括：
哨兵-3 温度与哨兵-2 NDVI 的散点图，作者图片
我们可以观察到此图中 NDVI 和温度之间的负相关关系，但某些点似乎是不规则值。这可能归因于图像中代表除植被和裸土之外特征的像素，例如开阔水域、城市区域，或者仅仅是质量较低的像素。为了去除这些像素，我们将只保留 NDVI 值在 0.1 到 0.6 之间，温度值在 300 开尔文到 330 开尔文之间，然后相应地更新图表：
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import linregress
# Filter the DataFrame based on the specified conditions
filtered_df = df_clean[(df_clean['NDVI'] >= 0.1) & (df_clean['NDVI'] <= 0.6) &
(df_clean['Temperature'] >= 300) & (df_clean['Temperature'] <= 330)]
# Create a scatter plot of Temperature versus NDVI
plt.figure(figsize=(8, 6))
sns.scatterplot(x='NDVI', y='Temperature', data=filtered_df, palette='coolwarm')
# Fit a linear regression model
slope, intercept = np.polyfit(filtered_df['NDVI'], filtered_df['Temperature'], 1)
# Plot the fitted line
x_line = np.linspace(min(filtered_df['NDVI']), max(filtered_df['NDVI']), 100)
y_line = slope * x_line + intercept
plt.plot(x_line, y_line, 'r', label='Fitted line')
# Create the equation string
equation_str = f'y = {slope:.2f}x + {intercept:.2f}'
# Display the equation on the plot
plt.text(min(filtered_df['NDVI']), max(filtered_df['Temperature']), equation_str, fontsize=12, color='red')
# Set plot title and axis labels
plt.title('1-1 Plot of Temperature vs NDVI', fontsize=16, fontweight='bold')
plt.xlabel('Temperature', fontsize=14)
plt.ylabel('NDVI', fontsize=14)
# Show the plot
plt.show()
输出结果将如下：
Sentinel-3 温度与 Sentinel-2 NDVI 的散点图，由作者绘制
现在，负相关关系更为明显，同时还有说明 NDVI 值与温度之间关系的方程。
📐 锐化热图像（1000 米到 10 米）
在这一步，我们将假设我们在 1000 米分辨率下发现的聚合 NDVI 和温度之间的方程也适用于 10 米分辨率。通过将方程应用于原始 NDVI 地图，我们可以估计 10 米分辨率的温度。
import rasterio
import numpy as np
import matplotlib.pyplot as plt
# Open the NIR and RED files
with rasterio.open('T10SFG_20230620T183919_B08_10m.jp2') as src:
nir = src.read(1)
meta = src.meta
with rasterio.open('T10SFG_20230620T183919_B04_10m.jp2') as src:
red = src.read(1)
# Calculate NDVI
ndvi = (nir - red) / (nir + red)
# Estimate temperature using NDVI
temp = -21.85 * ndvi + 314.9
# Create a color ramp for NDVI
ndvi_cmap = plt.cm.RdYlGn
# Plot NDVI and temperature side by side
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(12, 6))
# Plot NDVI
im1 = ax1.imshow(ndvi, cmap=ndvi_cmap, vmin=0, vmax=0.6)
ax1.set_title('NDVI', fontweight='bold', fontsize=14)
fig.colorbar(im1, ax=ax1, shrink=0.7)
# Plot temperature
im2 = ax2.imshow(temp, cmap=ndvi_cmap.reversed(), vmin=300, vmax=315)
ax2.set_title('Temperature', fontweight='bold', fontsize=14)
fig.colorbar(im2, ax=ax2, shrink=0.7)
plt.show()
在这里，我们读取近红外和红光波段，根据导出的方程计算（NDVI）并估计温度值。接下来，我们将 NDVI 和温度的结果以并排图的形式可视化。地图如下：
Sentinel-2 NDVI（10 米）位于左侧，下采样后的 Sentinel-3 温度（10 米）位于右侧。由作者可视化。
🗺️ 锐化热图像的可视化
在本节中，我们将进一步关注可视化方面。我们计划放大图像的中心区域，剪裁图像的中心部分，并展示 Sentinel-2 的 NDVI 地图、原始 Sentinel-3 热图像（1000 米分辨率）以及基于我们分析锐化的 Sentinel-3 图像（10 米）。这些可视化将并排展示以供比较和更详细的检查：
# Plot NDVI and temperature side by side
fig, axs = plt.subplots(ncols=3, figsize=(15, 5))
# Plot NDVI
vmin, vmax = 0, 0.6
ndvi_subset = ndvi[int(0.75 * ndvi.shape[0]):, int(0.75 * ndvi.shape[1]):]
im1 = axs[0].imshow(ndvi_subset, cmap=ndvi_cmap, vmin=vmin, vmax=vmax)
axs[0].set_title('NDVI', fontweight='bold', fontsize=14)
axs[0].set_xticks([])
axs[0].set_yticks([])
fig.colorbar(im1, ax=axs[0], shrink=0.7)
# Plot temperature
with rasterio.open('Sentinel-3_L2_LST_reproj_32610_clipped.tif') as src:
original_temp = src.read(1)
vmin, vmax = 300, 315
temp_subset = original_temp[int(0.75 * original_temp.shape[0]):, int(0.75 * original_temp.shape[1]):]
im3 = axs[1].imshow(temp_subset, cmap=ndvi_cmap.reversed(), vmin=vmin, vmax=vmax)
axs[1].set_title('Temperature', fontweight='bold', fontsize=14)
axs[1].set_xticks([])
axs[1].set_yticks([])
fig.colorbar(im3, ax=axs[1], shrink=0.7)
# Plot temperature
vmin, vmax = 300, 315
temp_subset = temp[int(0.75 * temp.shape[0]):, int(0.75 * temp.shape[1]):]
im2 = axs[2].imshow(temp_subset, cmap=ndvi_cmap.reversed(), vmin=vmin, vmax=vmax)
axs[2].set_title('Temperature', fontweight='bold', fontsize=14)
axs[2].set_xticks([])
axs[2].set_yticks([])
fig.colorbar(im2, ax=axs[2], shrink=0.7)
# Adjust the spacing between the subplots
fig.subplots_adjust(wspace=0.2)
plt.show()
地图如下：
Sentinel-2 NDVI（10 米）位于左侧，Sentinel-3 温度（1000 米）位于中间，下采样后的 Sentinel-3 温度（10 米）位于右侧，由作者可视化。
📄 结论
可见光和近红外（VNIR）波段与热图像之间的直接相关性已被证明是一种有用的方法，可以提高热图像的分辨率。这项技术在估计更高空间分辨率下的温度方面具有实际应用，尤其是在没有适当空间分辨率的卫星的情况下。这种方法在需要高分辨率热图并提供小尺度温度变化细节时是一个方便的工具。展望未来，更多配备先进热传感器的卫星的发射将使我们能够获得更高分辨率的更频繁热图像。在此之前，这种方法仍然是实现更高分辨率热图像的经济有效选项。
📚 参考文献
Copernicus Sentinel 数据[2024]用于 Sentinel 数据
Copernicus 服务信息[2024]用于 Copernicus 服务信息。
阿加姆，N.，库斯塔斯，W. P.，安德森，M. C.，李，F.，尼尔，C. M. U. (2007). 一种基于植被指数的热图像空间锐化技术。环境遥感，107(4)，545–558. ISSN 0034–4257.
高，F.，库斯塔斯，W. P.，安德森，M. C. (2012). 一种用于提高陆地热卫星图像清晰度的数据挖掘方法。遥感，4，3287–3319.
胡鲁纳，H.，科恩，Y.，卡尼利，A.，帕诺夫，N.，库斯塔斯，W. P.，阿加姆，N. (2019). 使用 Sentinel-2 视觉图像评估 TsHARP 在 Sentinel-3 卫星图像热增强中的效用。遥感，11，2304.
📱 在其他平台上与我联系，获取更多有趣的内容！领英，ResearchGate，GitHub，和***Twitter***。
这里是通过这些链接可以找到的相关帖子：
扫一扫
480
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
