深入挖掘ESA的Sentinel API
探索ESA Sentinel API中的深度
如何使用Python获取、分析和可视化卫星图像
本文中的所有图片均由作者创建。
欧洲空间局通过各种类型的遥感数据（如雷达和多光谱成像仪）支持欧洲联盟的地球观测组件Copernicus，其中包括Sentinel任务。
目前有六个正在进行中的Sentinel任务，其中三个可以通过它们的Python API轻松访问。以下是官方来源的引用：
Sentinel-1是一项极地轨道、全天候、日夜雷达成像任务，用于土地和海洋服务。Sentinel-1A于2014年4月3日发射，Sentinel-1B于2016年4月25日发射。这两个任务都是由Soyuz火箭从法属圭亚那的欧洲太空港发射进入轨道。Sentinel-1B任务于2022年结束，并计划尽快发射Sentinel-1C。
Sentinel-2是一项针对土地监测的极地轨道多光谱高分辨率成像任务，用于提供植被、土壤和水覆盖、内陆水道和沿海地区等成像。Sentinel-2还可以为紧急服务提供信息。Sentinel-2A于2015年6月23日发射，Sentinel-2B于2017年3月7日发射。
Sentinel-3是一项多仪器任务，用于精确可靠地测量海洋和陆地表面的海面地形、海洋和陆地表面温度、海洋和陆地颜色。该任务支持海洋预报系统以及环境和气候监测。Sentinel-3A于2016年2月16日发射，Sentinel-3B于2018年4月25日加入其孪生星进入轨道。
通过进一步了解，我们可以得知Sentinel-1的空间分辨率可达几米，而Sentinel-2的可视数据的最高分辨率为10米，而Sentinel-3的操作范围要大得多，根据传感器类型测量范围可达100公里。
好的，我们知道从哪里获取卫星数据了，而且看起来源（传感器）和空间分辨率的选择还非常丰富。可以指出，这只是冰山一角，因为这个卫星数据源列表还列出了其他相关资源。那么，我们使用这些不同类型的卫星数据做什么呢？首先，这里有50多个用例供您选择。
作为一个经验法则，我认为用例、问题的具体情况以及目标区域的地理空间特征和地形都是确定适合您的正确数据源的重要因素。然而，在日常实践中，据我的经验，以下几个主要因素起到关键作用：
- 价格（优先探索免费的Sentinel数据）
- 具有几米的空间分辨率，甚至可以捕捉到较小的城市结构
- 至少有几个波段，如可见光和近红外波段
- 时间频率
这些因素使得Sentinel-2可能是地理空间数据界使用最广泛的卫星数据源。在这些基础上，本文将向您展示如何获取Sentinel数据以及在下载时应该期望什么。我还将详细探讨图像记录和其中存储的信息的不同可能性和时间演变。
首先，按照官方文档和示例代码建立 API 连接。此外，我需要一个目标区域来下载图像。为了方便调试，我选择了我的家乡布达佩斯。我将使用 OSMNx 下载其行政边界。
import osmnx as ox # version: 1.0.1import matplotlib.pyplot as plt # version: 3.7.1city = 'Budapest'admin = ox.geocode_to_gdf(city)admin.plot()
现在点击 Sentinel API：
from sentinelsat import SentinelAPI, read_geojson, geojson_to_wkt # version 0.14# 要获取帐户，请在此注册：https://apihub.copernicus.eu/apihubuser = <添加您的用户名password = <添加您的密码>api = SentinelAPI(user, password, 'https://apihub.copernicus.eu/apihub') 要进行查询，最好有一个标明位置的平滑多边形。为此，我创建了布达佩斯的行政区域的凸包：
# 简化查询，提取输入多边形的凸包admin_polygon = admin.convex_hull.geometry.to_list()[0]admin_polygon输出：
在我们指定的位置、给定的时间范围和我们选择的航天器上搜索卫星图像。后者应为 Sentinel-A。此外，我们还可以根据云覆盖率进行过滤，这意味着如果图像太多云，我们可以立即放弃它。
# 这里我们可以指定位置（基于多边形）# 时间范围# 太空探测器# 和接受的云遮罩级别products = api.query(admin_polygon, date=('20150623', '20231006'), platformname='Sentinel-2', cloudcoverpercentage=(0, 100))len(products)正如这些单元格的输出所示，根据 Sentinel 文档，在 2015 年 06 月 23 日（任务开始）和 2023 年 10 月 06 日（我编写本文的时候），总共有 3876 张卫星图像记录与布达佩斯的行政边界重叠。我将云覆盖百分比设定为 100，表示没有基于云遮罩进行过滤。因此，我们应该有过去八年的所有图像标识符。
我还在这里说明，生成的产品列表包含所有卫星图像的标识符和元数据，但不包含图像本身。另外，如果我用 Sentinel-3 重复相同操作，会得到近 2 万条图像记录，但分辨率要低得多。
2. 探索元数据
让我们将产品列表转换为 Pandas DataFrame 并开始分析它！
import pandas as pd # version: 1.4.2products_gdf = api.to_geodataframe(products)products_gdf = products_gdf.sort_values(['beginposition'], ascending=[True])print(products_gdf.keys())print(len(products_gdf.keys()))products_gdf.head(3)此块的结果：
在计算由卫星图像标识符索引的表中的键数之后，人们可能已经对其中的41个特征列有了一个大致了解。
虽然这些字段中有很多技术信息，但我想更仔细地看一些特征。一方面，空间和时间维度通过生成日期、开始位置（作为日期时间信息）和几何形状（作为多边形、GIS数据类型）进行编码。另一方面，基于图像的一些有趣指标用于表征土地覆盖类型，包括我们在查询中已经看到的云覆盖百分比、植被百分比、水面百分比和雪冰百分比。这些环境指数是根据不同材料的不同光谱特性推导出的。注意：这些值都是捕捉整个瓷砖的总体平均值的汇总分数。关于这个主题的更多信息请参阅此处。
3. 空间维度
由于有几何维度，让我们看看它在地图上的样子！为此，我将通过可视化一组随机瓷砖来实现，经过几次运行后，这些瓷砖完全具有代表性。在可视化中，我使用了Folium和CartoDB Dark_Matter基础地图。
import foliumimport geopandas as gpdx, y = admin_polygon.centroid.xym = folium.Map(location=[y[0], x[0]], zoom_start=8, tiles='CartoDB Dark_Matter')# visualize a set of random tilespolygon_style = { 'fillColor': '#39FF14', 'color': 'black', 'weight': 3, 'opacity': 0}geojson_data = products_gdf[['geometry']].sample(10).to_json()folium.GeoJson( geojson_data, style_function=lambda feature: polygon_style).add_to(m)# add the admin boundaries on topadmin_style = {'fillColor': '#00FFFF', 'color': 'black','weight': 3, 'opacity': 100.0 }admin_geojson_data = admin[['geometry']].to_json()folium.GeoJson( admin_geojson_data, style_function=lambda feature: admin_style).add_to(m)# 显示地图m这段代码块的输出：
通过观察这个可视化图，很明显有一些瓷砖的区域出现重复。同时，可以看出，由于某种原因，有一些瓷砖将城市的行政边界分为两部分。这可能导致一个无法忍受的情况，即您希望完全分析您所需的目标区域的数据，但却发现它被分成两半。一个可能的解决方法是过滤掉与所需行政区域不完全重叠的瓷砖：
def compute_overlapping_area(tile, admin): return tile.intersection(admin_polygon).area / admin_polygon.areaproducts_gdf['overlapping_area_fraction'] = products_gdf.geometry.apply(lambda x: compute_overlapping_area(x, admin_polygon))products_gdf_f = products_gdf[products_gdf.overlapping_area_fraction==1]print(len(products_gdf))print(len(products_gdf_f))products_gdf_f.head(3)此单元格的结果：
通过应用此过滤器，我去除了大约一半的瓷砖。现在让我们看一看地图，看看它们与该城市的边界有多好重叠，没有任何瓷砖把城市分成两半：
import foliumimport geopandas as gpdx, y = admin_polygon.centroid.xym = folium.Map(location=[y[0], x[0]], zoom_start=8, tiles='CartoDB Dark_Matter')# visualize a set of random tilespolygon_style = { 'fillColor': '#39FF14', 'color': 'black', 'weight': 3, 'opacity': 0}geojson_data = products_gdf_f[['geometry']].sample(10).to_json()folium.GeoJson( geojson_data, style_function=lambda feature: polygon_style).add_to(m)# add the admin boundaries on topadmin_style = {'fillColor': '#00FFFF', 'color': 'black','weight': 3, 'opacity': 100.0 }admin_geojson_data = admin[['geometry']].to_json()folium.GeoJson( admin_geojson_data, style_function=lambda feature: admin_style).add_to(m)# 显示地图m
4. 数据集的时间维度
首先，让我们看看每天、每周和每月覆盖布达佩斯的图像数量。为了衡量时间，我将依赖 beginposition 字段。
# 假设 'beginposition' 是您的 GeoDataFrame 中的 Timestamp 列# 您可以将其转换为 DateTime 索引products_gdf_f_cntr = products_gdf_f.copy()products_gdf_f_cntr['beginposition'] = pd.to_datetime(products_gdf_f_cntr['beginposition'])products_gdf_f_cntr.set_index('beginposition', inplace=True)# 对数据进行重采样，统计每天、每周和每月的行数daily_counts = products_gdf_f_cntr.resample('D').count()weekly_counts = products_gdf_f_cntr.resample('W').count()monthly_counts = products_gdf_f_cntr.resample('M').count()fig, ax = plt.subplots(1, 3, figsize=(15, 5))for idx, (count_name, count_val) in enumerate([ ('每日计数', daily_counts), ('每周计数', weekly_counts), ('每月计数', monthly_counts), ]): ax[idx].plot(count_val.index[0:250], count_val['geometry'].to_list()[0:250]) ax[idx].set_xlabel('日期') ax[idx].set_ylabel('数量') ax[idx].set_title(count_name)plt.tight_layout()plt.suptitle('不同时间框架下拍摄的卫星图像数量', fontsize=20, y=1.15)plt.show()
这些图表展示了 Sentinel-2 探测器的前 250 天、前 250 周和前 250 个月（整个时间范围）。第一个图表显示每隔一天拍摄一次快照。第二个图表显示上一个图表的每周平均值，显示在最初的两年期间，卫星每周拍摄一到两次布达佩斯，然后从 2017 年到 2018 年，增加到每周 5-6 次。最后一个图表显示了整个时间范围内的相同趋势，以及在工作三年后，每月 25 张图像成为标准水平。
5. 土地覆盖变量的时间维度
现在，让我们来看一下植被百分比、水域百分比、雪冰百分比和云覆盖百分比的时间演变。正如前面的图像所示，早期的年份可能会显示不同的、可能是噪声的结果，因此我们要保持谨慎。在这里，我不会删除那些年份的数据，因为总共有八年的数据，删除其中的三年可能会丢失太多信息。首先，只是看一下每周聚合的原始值：
import pandas as pdimport matplotlib.pyplot as plt# 假设 'beginposition' 是您的 GeoDataFrame 中的 Timestamp 列# 您可以将其转换为 DateTime 索引products_gdf_f_cntr = products_gdf_f.copy()products_gdf_f_cntr['beginposition'] = pd.to_datetime(products_gdf_f_cntr['beginposition'])products_gdf_f_cntr.set_index('beginposition', inplace=True)# 对数据进行重采样，计算每周平均值weekly_averages = products_gdf_f_cntr[['vegetationpercentage', 'waterpercentage', 'snowicepercentage', 'cloudcoverpercentage']].resample('W').mean()# 创建包含四个子图的多图绘图fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, figsize=(10, 15))# 使用绿线绘制 'vegetationpercentage'ax1.plot(weekly_averages.index, weekly_averages['vegetationpercentage'], color='green', label='每周平均植被百分比')ax1.set_xlabel('日期')ax1.set_ylabel('百分比')ax1.set_title('每周平均植被百分比')ax1.legend()# 使用蓝线绘制 'waterpercentage'ax2.plot(weekly_averages.index, weekly_averages['waterpercentage'], color='blue', label='每周平均水域百分比')ax2.set_xlabel('日期')ax2.set_ylabel('百分比')ax2.set_title('每周平均水域百分比')ax2.legend()# 使用青色线绘制 'snowicepercentage'ax3.plot(weekly_averages.index, weekly_averages['snowicepercentage'], color='cyan', label='每周平均雪冰百分比')ax3.set_xlabel('日期')ax3.set_ylabel('百分比')ax3.set_title('每周平均雪冰百分比')ax3.legend()# 使用灰色线绘制 'cloudcoverpercentage'ax4.plot(weekly_averages.index, weekly_averages['cloudcoverpercentage'], color='gray', label='每周平均云覆盖百分比')ax4.set_xlabel('日期')ax4.set_ylabel('百分比')ax4.set_title('每周平均云覆盖百分比')ax4.legend()plt.tight_layout()plt.show()
每月聚合的情况如下：
这些时间序列告诉我们一些有趣的事情：
- 植被百分比清楚地显示季节性变化，每个春天一切都变绿，然后在秋天逐渐凋落，从50-60%下降到接近零。
- 相比之下，水域百分比在整个年份和观测期间都在0.8%左右波动。这是因为研究区域的地表水量非常少。尽管如此，水量似乎在冬季更频繁地下降，暗示一些淡水实体结冰。
- 至于雪，最明显的高峰发生在冬季，约为4-8%。尽管如此，根据我的个人经验，我可以告诉大家我们在大多数冬季中没有太多的雪。因此，仅测量1-2%的数值，尤其不是在冬季，可能会引起一些噪音甚至云的错误分类。
- 至于云，我们看到它们大多数时候与植被同步，遵循季节性模式。
这些措施之间的相关性在以下代码中也可见：
products_gdf_f_cntr[['vegetationpercentage', 'waterpercentage', 'snowicepercentage', 'cloudcoverpercentage']].corr()
6. 下载卫星图像
首先，针对Sentinel-2和Sentinel 3进行查询，选择今年8月的同一周，并尽量限制云覆盖。查看可用快照的数量：
# 查询瓷砖产品IDproducts_sent = api.query(admin_polygon, date=('20230806', '20230813'), platformname='Sentinel-2', cloudcoverpercentage=(0, 1))products_sent = api.to_geodataframe(products_sent)f, ax = plt.subplots(1,1,figsize=(6,4))admin.plot(ax=ax, color = 'none', edgecolor = 'k')ax.set_title('Sentinel-2，瓷砖数量 = ' + str(len(products_sent)))products_sent.plot(ax=ax, alpha = 0.3)# 过滤与布达佩斯不完全重叠的瓦片products_sent['overlapping_area_fraction'] = products_sent.geometry.apply(lambda x: compute_overlapping_area(x, admin_polygon))products_sent = products_sent[products_sent.overlapping_area_fraction==1]f, ax = plt.subplots(1,1,figsize=(6,4))admin.plot(ax=ax, color = 'none', edgecolor = 'k')ax.set_title('Sentinel-2，瓷砖数量 = ' + str(len(products_sent)))products_sent.plot(ax=ax, alpha = 0.3)len(products_sent)此代码块的结果：
现在根据图像数据下载
# 下载第一批作为卫星图像product_ids = products_sent.index.to_list()for prod in product_ids: api.download(prod)注意 — 在这种情况下，您可能会收到此通知，请等几个小时然后再运行下载程序。
产品a3c61497-d77d-48da-9a4d-394986d2fe1d不在线。触发从长期存档检索。7。打开和可视化已下载的图像
在这里，您将找到有关数据格式的详细描述，并有关于文件夹结构的非常好看的可视化。一旦打开图像目录，您可能会找到不同的波段。每个波段的含义以及其空间分辨率在这篇文章中都有很好的总结，因为13个波段的空间分辨率范围是10到60米。一些亮点：
- 蓝色（B2），绿色（B3），红色（B4）和近红外或NIR（B8）通道具有10米的分辨率。
- 然后，接下来，它的植被红边（B5），近红外NIR（B6，B7和B8A）以及短波红外SWIR（B11和B12）具有10米的分辨率。
- 最后，其沿海气溶胶（B1）和SWIR卷云波段（B10）具有60米的像素大小。
这就是它的样子：
#在解压下载的文件夹后：import osimage_path = 'S2B_MSIL1C_20230810T094549_N0509_R079_T34TCT_20230810T124346.SAFE/GRANULE/L1C_T34TCT_A033567_20230810T095651/IMG_DATA'sorted(os.listdir(image_path))
当使用rasterio可视化B4，即红色波段时，整个图块的外观如下：
import rasteriofrom rasterio.plot import showimage_file = 'T34TCT_20230810T094549_B04.jp2' with rasterio.open(image_path + '/' + image_file) as src: image = src.read(1) # 需要时更改波段索引 plt.figure(figsize=(10, 10)) plt.imshow(image, cmap='Reds') # 您可以更改颜色映射 plt.title(image_file) plt.colorbar() plt.show()输出：
现在将重点放在布达佩斯，并将整个栅格图块按城市的行政边界遮罩，逐个遮罩处理R、G和B波段：
from rasterio import maskf, ax = plt.subplots(1,3,figsize=(15,5))for idx, (band_name, band_num, color_map) in enumerate([('蓝色', 'B02', '蓝色'), ('绿色', 'B03', '绿色'), ('红色', 'B04', '红色')]): raster_path = image_path + '/T34TCT_20230810T094549_' + band_num + '.jp2' with rasterio.open(raster_path) as src: polygons = admin.copy().to_crs(src.crs) geom = polygons.geometry.iloc[0] masked_image, _ = mask.mask(src, [geom], crop=True) ax[idx].imshow(masked_image[0], cmap=color_map) ax[idx].set_title('布达佩斯Sentinel 2 - ' + band_name + '波段')
最后，让我们将它们组合起来，获得一张布达佩斯的RGB图像。首先，我将完整的切片组合成一张RGB图像，然后读取它，按照官方说明进行直方图均衡化，最后得到最终的图。
# 获取波段位置band_blue = '/T34TCT_20230810T094549_B02.jp2'band_green = '/T34TCT_20230810T094549_B03.jp2'band_red = '/T34TCT_20230810T094549_B04.jp2'# 读取波段并创建完整的RGB切片b2 = rasterio.open(image_path + '/' + band_blue)b3 = rasterio.open(image_path + '/' + band_green)b4 = rasterio.open(image_path + '/' + band_red)# 将完整切片导出为tif文件meta = b4.metameta.update({"count": 3})prefire_rgb_path = 'budapest_rgb.tif'with rasterio.open(prefire_rgb_path, 'w', **meta) as dest:dest.write(b2.read(1),1)dest.write(b3.read(1),2)dest.write(b4.read(1),3)# 读取并显示裁剪版本from skimage import exposureimg = rasterio.open('budapest_rgb_cropped.tif')image = np.array([img.read(3), img.read(2), img.read(1)])image = image.transpose(1,2,0)# 进行直方图均衡化p2, p98 = np.percentile(image, (2,98))image = exposure.rescale_intensity(image, in_range=(p2, p98)) / 100000f, ax = plt.subplots(1,1,figsize=(15,15))rasterio.plot.show(image.transpose(2,0,1), transform=img.transform, ax = ax)ax.axis('off')plt.savefig('budapest_rgb_cropped.png', dpi = 700, bbox_inches = 'tight')输出：
结束语
快速总结一下，看看本文发生了什么：
- 关于Sentinel卫星平台的简要概述
- 如何查询多个图像标识符，包括它们的元数据的示例
- 如何简单地基于存储和定向在元数据中的切片的聚合信息进行时间分析
- 如何下载、存储和可视化单个图像
所有这些步骤旨在将卫星图像处理和分析添加到您日常的地理空间数据科学工具堆栈中，它可以覆盖从城市规划到环境监测和农业的许多用例。
