来源：ChatGPT 4-o
所有代码均使用 Python 编写，并在 Google Colab 中进行测试。
目录
-
🌟 简介
-
📌 研究区域 (AOI)
-
💾 加载 Sentinel-2 影像
-
⏳ 从卫星影像中提取时间序列
-
🌍 开发具有时间序列的交互式地图
-
📄 结论
-
📚 参考文献
🌟 简介
一段时间以来，我一直寻找一种简单直接的方法来创建一个交互式地图，当用户点击特定位置时显示时间序列数据。我探索了 folium 库，并找到了如何将卫星影像提取的时间序列数据库与 folium 连接起来实现这一功能的方法。今天，我将分享我的方法。
在这篇文章中，我将编写两个函数。第一个函数在不下载的情况下加载卫星数据，第二个函数提取数据和时间戳，创建一个数据格式的时间序列。然后，我们将遍历我们的研究区域（AOI），运行这两个函数并提取我们的 AOI 的时间序列数据。最后，我们将使用生成的时间序列数据和 folium 库在交互式地图上地理空间地显示它。
在本篇文章结束时，您将能够将提取的任何变量或参数的时间序列数据与交互式地图连接起来，并在地理空间上可视化它们。作为一个例子，我将提取加利福尼亚湖泊表面积的时间序列并在交互式地图上显示它们。然而，您可以根据这些步骤对任何您感兴趣的变量进行操作，例如多云天数、云层覆盖、洪水区域、植被指数、降水量、空气温度等。如果您对此感兴趣或一直在寻找这些技巧和窍门，请继续阅读！
📌 研究区域 (AOI)
如简介中所述，您可以根据以下步骤在交互式地图上显示任何变量在任何位置的时间序列数据。对于这个例子，我将计算加利福尼亚几个湖泊的水像素数量，并使用 2024 年捕获的所有 Sentinel-2 影像计算它们的表面积。我在 QGIS 中围绕这些湖泊绘制了多边形，并将它们保存为 shapefile。如果您想学习如何为研究区域（AOI）创建 shapefile，请参考我在 Medium 上的第一篇文章中标题为**“🛠️ 创建 QGIS 中的 Shapefile”**的部分：
下面是我在 QGIS 中围绕湖泊绘制的多边形的快照：
QGIS 中湖泊周围的多边形，图像由作者提供
💾 加载 Sentinel-2 图像
本节的目标是在不下载的情况下将存档卫星图像加载到内存中。下载特定区域长时间段的卫星图像既耗时，又计算量大，效率低下，尤其是如果你想从小区域中探索整个场景的变化。
为了克服这些问题，我写了一篇帖子，展示了如何仅用 12 行代码加载卫星图像而不下载它们，你可以在这里找到：
在本节中，我们将使用该帖子中展示的相同模板编写一个函数。使用此函数，你可以轻松加载特定时间段内覆盖你的 AOI 的卫星图像，无论这个时间段是长是短：
from pystac_client import Client
from odc.stac import load
def search_satellite_images(collection="sentinel-2-l2a",
bbox=[-120.15,38.93,-119.88,39.25],
date="2023-01-01/2023-03-12",
cloud_cover=(0, 10)):
"""
Search for satellite images based on collection, bounding box, date range, and cloud cover.
:param collection: Collection name (default: "sentinel-2-l2a").
:param bbox: Bounding box [min_lon, min_lat, max_lon, max_lat] (default: Lake Tahoe Region).
:param date: Date range "YYYY-MM-DD/YYYY-MM-DD" (default: "2023-01-01/2023-12-30").
:param cloud_cover: Tuple representing cloud cover range (min, max) (default: (0, 10)).
:return: Data loaded based on search criteria.
"""
# Define the search client
client=Client.open("https://earth-search.aws.element84.com/v1")
search = client.search(collections=[collection],
bbox=bbox,
datetime=date,
query=[f"eo:cloud_cover<{cloud_cover[1]}", f"eo:cloud_cover>{cloud_cover[0]}"])
# Print the number of matched items
print(f"Number of images found: {search.matched()}")
data = load(search.items(), bbox=bbox, groupby="solar_day", chunks={})
print(f"Number of days in data: {len(data.time)}")
return data
此函数允许你指定搜索卫星图像的参数，使其灵活且易于使用，适用于不同的 AOI 和时间段。它根据集合、边界框、日期范围和云层覆盖等标准搜索卫星图像。它使用 Earth Search API 查找图像，并打印匹配的数量，以立方体格式返回裁剪图像。
⏳ 从卫星图像中提取时间序列
现在我们已经有了用于加载我们的 AOI 裁剪图像的函数，我们需要定义第二个函数来提取我们所需的信息，包括图像，并将其放入 DataFrame 中，以便作为时间序列数据库在下一步中使用（在地图上显示）。再次强调，你可以提取你需要的数据，但我认为看看 2024 年 Sentinel-2 捕获的所有图像中主要加利福尼亚湖泊的表面积会很有趣，包括非常近期的图像。
要实现这一点，我们可以从 Sentinel-2 图像的场景分类层中提取被分类为水像素的像素。换句话说，我们需要计算每个场景中的水像素数量。鉴于像素分辨率为 10m x 10m，将计数乘以 100 平方米（10m x 10m）将给出每个湖泊的表面积。然而，在这个特定例子中，我们需要确保卫星捕获的图像覆盖了每个场景中的整个湖泊。为了说明这一点，让我们看看这些图像，它们是在 1 月 7 日和 1 月 4 日捕获的，覆盖了其中一个湖泊。
作者可视化于 1 月 4 日在清水湖上捕获的 Sentinel 图像
Sentinel 图像于 1 月 7 日在 Clear Lake 上捕获，由作者可视化
如您所见，在一幅图像中，整个湖泊都包含在场景中，而在另一幅中，只有 10%的湖泊被覆盖。如果我们不考虑这一点来计算水像素，我们可能会错误地得出湖泊表面积显著下降的结论，而实际上，这只是由于覆盖不完整。为了解决这个问题，我们需要计算我所说的“覆盖比”。在代码中，我跳过了覆盖比低于 80%的图像的分析。您可以随意将其更改为 100%。
因此，第二个函数将是：
import pandas as pd
import numpy as np
def count_water_pixels(data,lake_id):
"""
Counts water pixels from Sentinel-2 SCL data for each time step.
:param data: xarray Dataset with Sentinel-2 SCL data.
:return: DataFrame with dates, water counts, and snow counts.
"""
water_counts = []
date_labels = []
water_area =[]
coverage_ratio=[]
# Determine the number of time steps
numb_days = len(data.time)
# Iterate through each time step
for t in range(numb_days):
scl_image = data[["scl"]].isel(time=t).to_array()
dt = pd.to_datetime(scl_image.time.values)
year = dt.year
month = dt.month
day = dt.day
date_string = f"{year}-{month:02d}-{day:02d}"
print(date_string)
# Count the number of pixels corresponding to water
count_water = np.count_nonzero(scl_image == 6) # Water
surface_area = count_water*10*10/(10**6)
count_pixels = np.count_nonzero((scl_image == 1) | (scl_image == 2) | (scl_image == 3)|(scl_image == 4) | (scl_image == 5) |
(scl_image == 6) | (scl_image == 7) |(scl_image == 8) |(scl_image == 9) |(scl_image == 10) |(scl_image == 11))
total_pixels = data.dims['y']*data.dims['x']
coverage = count_pixels*10*10/1e6
lake_area = total_pixels*10*10/1e6
ratio=coverage/lake_area
print(coverage)
print(lake_area)
print(ratio)
if ratio <0.8:
continue
# Append
water_counts.append(count_water)
date_labels.append(date_string)
water_area.append(surface_area)
coverage_ratio.append(ratio)
# Convert date labels to pandas datetime format
datetime_index = pd.to_datetime(date_labels)
# Create a dictionary for constructing the DataFrame
data_dict = {
'Date': datetime_index,
'ID': lake_id,
'Water Counts': water_counts,
'Pixel Counts': count_pixels,
'Total Pixels': total_pixels,
'Coverage Ratio': coverage_ratio,
'Water Surface Area': water_area
}
# Create the DataFrame
df = pd.DataFrame(data_dict)
return df
此函数遍历数据集中的每个时间步，计算水像素，计算表面积，并计算覆盖比。如果覆盖比低于 80%，则跳过该时间步。然后，该函数将计数、日期、表面积和覆盖比追加到列表中，并返回一个包含这些值的 DataFrame，以及湖泊 ID 和总像素计数。
如果您想了解更多关于覆盖问题及其解决技巧的信息，请参阅我在这里写的文章的（📈 从统计文件中获取大盐湖面积的时间序列）部分：
🌍 开发具有时间序列的交互式地图
在本节中，我们将编写三个脚本。第一个脚本是一个函数，用于提取我们多边形（AOI）的边界框以及中心坐标。我们需要边界框来运行我们的第一个函数（search_satellite_images），以及中心坐标来在我们地图上标记湖泊。这可以通过以下代码完成：
import geopandas as gpd
import pandas as pd
def get_centroids_and_bboxes(shapefile_path):
"""
Processes a shapefile to return a DataFrame containing the ID, centroid,
and bounding box (bbox) of each polygon.
:param shapefile_path: Path to the shapefile.
:return: A pandas DataFrame with the ID, centroid, and bbox of each polygon.
"""
# Load the shapefile
gdf = gpd.read_file(shapefile_path)
# Reproject to EPSG:4326
gdf_proj = gdf.to_crs("EPSG:4326")
centroids = []
bboxes = []
# Process each polygon to get its centroid and bbox
for index, row in gdf_proj.iterrows():
# Centroid
centroid_lat = row.geometry.centroid.y
centroid_lon = row.geometry.centroid.x
centroids.append((centroid_lat, centroid_lon))
# Bounding Box
minx, miny, maxx, maxy = row.geometry.bounds
bbox = (minx, miny, maxx, maxy)
bboxes.append(bbox)
# Create the DataFrame
df = pd.DataFrame({
'ID': gdf_proj.index,
'Centroid_Lat': [lat for lat, lon in centroids],
'Centroid_Lon': [lon for lat, lon in centroids],
'BBox_Min_Lon': [bbox[0] for bbox in bboxes],
'BBox_Min_Lat': [bbox[1] for bbox in bboxes],
'BBox_Max_Lon': [bbox[2] for bbox in bboxes],
'BBox_Max_Lat': [bbox[3] for bbox in bboxes]
})
return df
shapefile_path = "lakes_boundry.shp"
lakes_df = get_centroids_and_bboxes(shapefile_path)
print(lakes_df)
如果您遵循步骤并成功运行代码，您应该会看到以下格式的类似 DataFrame，用于您的多边形：
AOI 边界框和中点坐标的 DataFrame，图像由作者提供
第二个脚本涉及遍历 Sentinel-2 在 2024 年捕获的所有湖泊图像，如果覆盖比高于 80%，则运行第二个函数来计算每幅图像的表面积，并报告一个显示每个湖泊表面积的 DataFrame 时间序列：
import matplotlib.pyplot as plt
import pandas as pd
from datetime import datetime
import numpy as np
all_water_pixels_dfs = []
for lake_id in lakes_df.ID:
print(lake_id)
lake_df = lakes_df[lakes_df['ID'] == lake_id]
if not lake_df.empty:
bbox = [lake_df.iloc[0].BBox_Min_Lon, lake_df.iloc[0].BBox_Min_Lat,
lake_df.iloc[0].BBox_Max_Lon, lake_df.iloc[0].BBox_Max_Lat]
data = search_satellite_images(collection="sentinel-2-l2a",
date="2024-01-01/2024-05-14",
cloud_cover=(0, 5),
bbox=bbox)
# Pass the lake_id
water_pixels_df = count_water_pixels(data, lake_id)
# Append
all_water_pixels_dfs.append(water_pixels_df)
# Concatenate all DataFrames into a single DataFrame
final_df = pd.concat(all_water_pixels_dfs, ignore_index=True)
最终的 DataFrame 将如下所示，总结图像日期、水像素数量、总像素数、覆盖比和表面积：
基于 Sentinel-2 在 2024 年的图像显示每个湖泊表面积时间序列的最终 DataFrame，图像由作者提供
我们几乎就完成了！
只需再走一步，就能在地图上看到时间序列。现在我们有了表面面积的时间序列数据，我们可以使用 Folium 库展示两个东西：(1) 湖泊在地图上的中心点，以及(2) 当用户点击每个湖泊时弹出的表面面积时间序列图表。这可以通过以下代码实现：
import folium
import plotly.express as px
import os
# Function to plot time series and save as HTML
def plot_timeseries_for_spot(spot_id, ts_df):
df_spot = ts_df[ts_df['ID'] == spot_id]
print(df_spot)
fig = px.line(df_spot, x='Date', y='Water Surface Area', title=f'Time Series for Lake {spot_id}')
# Add X and Y axis labels
fig.update_layout(
xaxis_title="Date",
yaxis_title="Water Surface Area (sq km)"
)
filepath = f'tmp_{int(spot_id)}.html'
fig.write_html(filepath, include_plotlyjs='cdn')
return filepath
# Create a map
m = folium.Map(location=[35.5, -119.5], zoom_start=7)
# Add markers with Plotly time series popups
for index, row in lakes_df.iterrows():
html_path = plot_timeseries_for_spot(row['ID'], final_df)
iframe = folium.IFrame(html=open(html_path).read(), width=500, height=300)
popup = folium.Popup(iframe, max_width=2650)
folium.Marker([row['Centroid_Lat'], row['Centroid_Lon']], popup=popup).add_to(m)
m.save('map_with_timeseries.html')
# Clean up temporary HTML files
for spot_id in lakes_df['ID']:
os.remove(f'tmp_{spot_id}.html')
在那个脚本中，该函数过滤每个湖泊的时间序列数据，使用 Plotly 生成线形图，并将图表保存为 HTML 文件。然后，它使用 Folium 初始化地图。接着，它遍历湖泊 DataFrame，为每个湖泊的中心坐标添加标记，并为每个标记附加一个弹出窗口，显示时间序列图表。最终的地图保存为 HTML 文件。最后，它删除为 Plotly 图表创建的临时 HTML 文件以进行清理。
完成！
如果您打开内容文件夹中创建的 HTML 文件，您将看到每个湖泊在地图上标记的中心坐标。
在地图上标记了 AOI（本例中的湖泊）的中心，图片由作者提供
让我们点击每个湖泊，看看时间序列是否会弹出：
2024 年从 Sentinel-2 图像中计算出的湖#5 的表面面积时间序列，图片由作者提供。
2024 年从 Sentinel-2 图像中计算出的湖#7 的表面面积时间序列，图片由作者提供。
经过这些努力后，终于得到了这样一张实用的地图，不是吗？😀
📄 结论
几乎每个月我都会看到新的包和库被设计出来，用于以实用方式提取、分析、显示和可视化数据。然而，在这个领域仍然存在两个挑战。第一个是正确分析数据需要足够的经验，以确保从千兆或太字节的数据中提取的数据是准确的。第二个是创建一个架构，将这些库连接起来，以产生有意义且实用的东西。
在处理图像时，我试图强调简单的错误如何导致处理数据中的重大错误。在可视化部分，我展示了如何连接 Folium、Plotly 和用于提取裁剪图像的新 API，可以创建一个用于使用遥感观测监测各种现象的有用工具。希望您喜欢阅读，如果您有任何问题，请随时联系。
📚 参考文献
github.com/stac-utils/pystac-client/blob/main/docs/quickstart.rst
www.element84.com/earth-search/examples/
Copernicus Sentinel 数据[2024]用于 Sentinel 数据
Copernicus 服务信息[2024]用于 Copernicus 服务信息。
📱 在其他平台上与我联系，获取更多引人入胜的内容！领英，ResearchGate，GitHub，以及***Twitter***。
这里是通过以下链接可用的相关帖子：
一个长 Python 脚本，用于使用地理空间大数据制作短视频
扫一扫
1317
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
