创建一个交互式地图以显示卫星影像的时间序列

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

towardsdatascience.com/create-an-interactive-map-to-display-time-series-of-satellite-imagery-e9346e165e27

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/4e65eddcb6e836c919ed8b36c508ff1b.png

来源：ChatGPT 4-o

所有代码均使用 Python 编写，并在 Google Colab 中进行测试。

简介

长期以来，我一直在寻找一种简洁高效的方式，来构建一个交互式地图——当用户点击地图上的某个位置时，能够动态显示该地点的时间序列数据。经过对 folium 库的深入探索，我成功实现了将从卫星影像中提取的时间序列信息与地图进行集成的方法。本文将详细介绍这一过程。

我们将定义两个核心函数：第一个用于直接从云端加载 Sentinel-2 卫星影像而无需本地下载；第二个则负责提取影像中的时间戳和对应数据，整理成标准的时间序列格式。随后，程序会遍历研究区域（AOI）内的各个目标点或区域，调用这两个函数获取每个位置的时间序列数据。最终，利用 folium 将这些时间序列结果以地理空间可视化的方式呈现在交互式地图上。

通过本文的学习，你将掌握如何将任意变量（如植被指数、降水、地表温度等）的时间序列数据与地图联动展示。作为示例，本文将以加利福尼亚州多个湖泊为例，计算其表面积随时间的变化趋势并实现点击查询功能。但该方法可广泛应用于其他遥感监测场景，包括洪水范围变化、云频统计、干旱评估等。

研究区域 (AOI)

正如前文所述，该流程适用于任何地理位置与感兴趣参数的组合。本案例聚焦于美国加利福尼亚州境内的若干湖泊，目标是基于 2024 年全年获取的 Sentinel-2 影像，估算各湖泊水体像素数量，进而推算其表面积变化趋势。

为精确定义分析范围，我在 QGIS 软件中手动绘制了围绕这些湖泊的多边形边界，并将其保存为 shapefile 格式文件，以便后续在代码中调用。以下为绘制完成的多边形示意图：

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/ae52dfa15963bf084a122203ba75b515.png

QGIS 中湖泊周围的多边形，图像由作者提供

加载 Sentinel-2 影像

本阶段的核心目标是在不执行本地下载的前提下，将远端存档的卫星影像直接加载至运行内存中。传统方式中，若需获取某区域长时间跨度的遥感数据，往往涉及大量影像的下载、存储与处理，不仅耗时且资源消耗巨大，尤其对于小范围 AOI 来说效率极低。

为此，我此前已发布一篇技术文章，介绍了一种仅用 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)):
    """
    基于指定的数据集、地理范围、时间区间和云覆盖率筛选卫星影像。
    
    :param collection: 数据集名称（默认："sentinel-2-l2a"）
    :param bbox: 地理边界框 [最小经度, 最小纬度, 最大经度, 最大纬度]（示例为太浩湖区域）
    :param date: 时间范围，格式为 "YYYY-MM-DD/YYYY-MM-DD"
    :param cloud_cover: 允许的云覆盖率区间，取值范围 0–100 的元组
    """

从卫星影像中提取时间序列

在成功加载影像后，下一步是从每一景影像中提取有效观测值及其对应的时间标签，形成结构化的时间序列数据。这一步通常包括影像预处理、波段计算（如归一化水体指数 NDWI）、掩膜应用（去除云、阴影、非水体干扰）以及按 AOI 进行区域统计。

通过对每一张符合条件的影像重复上述操作，我们便可获得按时间排序的序列化输出，例如每个日期对应的湖泊水体像元数、平均反射率、面积估算值等。

开发具有时间序列的交互式地图

完成时间序列数据的提取后，接下来的关键步骤是将其与交互式地图集成。我们采用 Python 的 folium 库构建底图，并在每个 AOI 位置添加标记点。当用户点击某一标记时，前端将触发回调事件，动态渲染该位置对应的时间序列图表（如折线图或柱状图），实现“点击即见趋势”的交互体验。

这种结合地理空间定位与时间维度分析的可视化方式，极大提升了遥感数据的可解释性与实用性，特别适合环境监测、城市规划、农业管理等领域。

结论

本文介绍了一套完整的流程，用于实现遥感时间序列数据与交互式地图的无缝集成。通过云端加载 Sentinel-2 影像、自动化提取 AOI 内的时间序列信息，并结合 folium 实现地理可视化，用户可以轻松监控任意区域、任意变量的动态变化。

该方法避免了大规模数据下载带来的开销，提高了处理效率，同时具备良好的扩展性和复用性。无论是湖泊萎缩、植被生长周期还是城市扩张分析，均可依此模式快速搭建分析系统。

参考文献

（原文未提供具体参考文献列表，此处保留标题结构）

该函数用于定义搜索卫星影像的参数，具有较高的灵活性和易用性，适用于不同区域（AOI）及时间范围的数据获取。通过指定数据集合、地理边界框（bbox）、时间范围以及云覆盖率范围等条件，利用 Earth Search API 进行影像检索。程序会输出符合条件的影像数量，并以立方体格式返回裁剪后的数据。

函数内部实现如下：

:param cloud_cover: 表示云覆盖率范围的元组 (最小值, 最大值)（默认为 (0, 10)）
:return: 根据搜索条件加载的数据
"""
# 初始化搜索客户端
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(f"找到的影像数量: {search.matched()}")
data = load(search.items(), bbox=bbox, groupby="solar_day", chunks={})
print(f"数据中包含的天数: {len(data.time)}")
return data

towardsdatascience.com/create-an-interactive-map-to-display-time-series-of-satellite-imagery-e9346e165e27

从卫星影像中提取时间序列数据

在实现了对目标区域（AOI）影像加载的功能后，下一步是构建一个辅助函数，用于从这些影像中提取关键信息，并将其组织成可用于后续分析的 DataFrame 结构，例如作为时间序列数据用于地图展示或其他统计用途。

虽然用户可以根据需求自定义提取内容，但本文将以美国加利福尼亚州主要湖泊为例，分析2024年期间Sentinel-2卫星所捕获的所有图像中水体表面积的变化趋势，包括最新获取的影像数据。

具体方法是利用Sentinel-2影像自带的场景分类层（Scene Classification Layer, SCL），识别出被标记为“水体”的像素点。由于Sentinel-2的分辨率为10米×10米，每个像素对应100平方米的地表面积，因此将水体像素总数乘以100即可估算出湖泊的表面积。

github.com/stac-utils/pystac-client/blob/main/docs/quickstart.rst

然而，需要注意的是，并非所有影像都完整覆盖了整个湖泊。如以下两幅图像所示，分别拍摄于1月4日和1月7日，对同一湖泊的覆盖程度存在显著差异：

• 一幅图像完整包含了整个湖泊；
• 另一幅则仅覆盖了约10%的湖面区域。

如果不考虑这种空间覆盖的不一致性，直接计算水体像素数量，可能会误判为湖泊面积急剧缩小，而实际上只是影像覆盖范围不同所致。为此，引入“覆盖比”这一指标进行筛选。在实际编码过程中，仅保留覆盖比高于80%的影像参与分析（也可根据需要调整至更高阈值，如100%）。

基于上述逻辑，定义第二个核心函数如下：

import pandas as pd
import numpy as np

def count_water_pixels(data, lake_id):
    """
    针对每一步时间数据，统计Sentinel-2 SCL中的水体像素数量。
    :param data: 包含Sentinel-2 SCL数据的xarray Dataset
    :return: 包含日期、水体像素计数和积雪像素计数的DataFrame
    """
    water_counts = []
    date_labels = []
    water_area = []
    coverage_ratio = []

    # 获取总的时间步长数量
    numb_days = len(data.time)

该函数通过遍历数据集中的每一个时间步，执行以下操作：统计水体像素数量、计算水体表面积以及影像覆盖比例。当覆盖比例低于80%时，跳过当前时间步的数据；否则，将水体像素数、对应日期、水体面积和覆盖比等信息依次添加到相应的列表中。

在循环结束后，将日期标签转换为pandas的datetime格式，并构建一个包含如下字段的数据字典：

• Date：观测日期
• ID：湖泊标识符（lake_id）
• Water Counts：检测到的水体像素数量
• Pixel Counts：有效像素总数
• Total Pixels：区域总像素数
• Coverage Ratio：有效影像覆盖比例
• Water Surface Area：水体表面面积（单位：平方公里）

随后，利用该字典创建一个pandas DataFrame并返回结果。此过程实现了对湖泊水域随时间变化的量化分析。

towardsdatascience.com/create-an-interactive-map-to-display-time-series-of-satellite-imagery-e9346e165e27

从统计文件提取大盐湖面积时间序列

为了获取大盐湖的历史面积变化情况，我们可基于处理后的卫星影像结果生成连续的时间序列数据。该时间序列可用于分析湖泊萎缩趋势，支持生态环境监测与水资源管理。

开发带时间序列的交互式地图

接下来我们将编写三个独立脚本以实现可视化功能。第一个脚本用于提取研究区域（AOI）多边形的边界框及其中心坐标。其中，边界框用于调用影像搜索函数 search_satellite_images，而中心坐标则用于在地图上标注湖泊位置。

以下是实现该功能的代码：

import geopandas as gpd
import pandas as pd

def get_centroids_and_bboxes(shapefile_path):
    """
    读取矢量文件并返回每个多边形的ID、几何中心点和边界框。
    
    :param shapefile_path: 矢量文件路径
    :return: 包含ID、centroid和bbox的DataFrame
    """
    # 读取shapefile
    gdf = gpd.read_file(shapefile_path)
    
    # 计算中心点（经纬度）
    gdf['centroid'] = gdf.geometry.centroid
    gdf['lon'] = gdf.centroid.x
    gdf['lat'] = gdf.centroid.y
    
    # 提取边界框 [minx, miny, maxx, maxy]
    gdf['bbox'] = gdf.geometry.bounds.apply(lambda row: 
              [row['minx'], row['miny'], row['maxx'], row['maxy']], axis=1)
    
    # 构建输出数据框
    result_df = gdf[[gdf.columns[0], 'lon', 'lat', 'bbox']].copy()
    result_df.columns = ['ID', 'lon', 'lat', 'bbox']
    
    return result_df

该函数返回一个包含湖泊ID、经度、纬度和空间范围（bbox）的信息表，为后续影像检索与地图标注提供基础支持。

def get_centroids_and_bboxes(shapefile_path):
    """
    读取Shapefile文件，计算每个面状要素的中心点与边界框，并返回包含这些信息的数据框。
    :param shapefile_path: Shapefile 文件路径
    :return: 包含 ID、中心点坐标和边界框的 pandas DataFrame
    """
    # 加载地理矢量数据
    gdf = gpd.read_file(shapefile_path)
    
    # 将坐标系转换为 WGS84（EPSG:4326）
    gdf_proj = gdf.to_crs("EPSG:4326")
    
    centroids = []
    bboxes = []

    # 遍历每个多边形要素，提取其中心点和包围盒
    for index, row in gdf_proj.iterrows():
        # 获取中心点经纬度
        centroid_lat = row.geometry.centroid.y
        centroid_lon = row.geometry.centroid.x
        centroids.append((centroid_lat, centroid_lon))

        # 提取几何对象的最小外接矩形（minx, miny, maxx, maxy）
        minx, miny, maxx, maxy = row.geometry.bounds
        bbox = (minx, miny, maxx, maxy)
        bboxes.append(bbox)

    # 构建结果数据表
    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 路径
shapefile_path = "lakes_boundry.shp"
lakes_df = get_centroids_and_bboxes(shapefile_path)
print(lakes_df)

若代码正确执行并按步骤运行，您将获得一个类似如下结构的 DataFrame，展示各多边形区域的 AOI 边界框及其中心点坐标：

towardsdatascience.com/create-an-interactive-map-to-display-time-series-of-satellite-imagery-e9346e165e27

import matplotlib.pyplot as plt
import pandas as pd
from datetime import datetime
import numpy as np

# 存储所有湖泊水体像元统计结果的列表
all_water_pixels_dfs = []

# 遍历每一个湖泊 ID
for lake_id in lakes_df.ID:
    print(lake_id)
    
    # 根据 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
        ]
        
        # 查询指定时间段内符合条件的 Sentinel-2 图像（云量 0-5%）
        data = search_satellite_images(
            collection="sentinel-2-l2a",
            date="2024-01-01/2024-05-14",
            cloud_cover=(0, 5),
            bbox=bbox
        )
        
        # 计算该湖泊在每景影像中的水体像素数量，传入 lake_id 标识
        water_pixels_df = count_water_pixels(data, lake_id)
        
        # 将单个湖泊的结果加入总列表
        all_water_pixels_dfs.append(water_pixels_df)

# 合并所有湖泊的时间序列数据为单一 DataFrame
final_df = pd.concat(all_water_pixels_dfs, ignore_index=True)

最终生成的 DataFrame 将包含图像日期、水体像素数量、总像素数、覆盖比例以及湖泊表面积等信息，如下所示：

towardsdatascience.com/create-an-interactive-map-to-display-time-series-of-satellite-imagery-e9346e165e27

该数据基于 2024 年 Sentinel-2 卫星影像构建，展示了各个湖泊表面积随时间变化的时间序列。DataFrame 的内容由作者整理提供。

接下来只剩下最后一步：将时间序列可视化在地图上。目前我们已经获得了湖泊表面积的时序数据，可以借助 Folium 库实现以下两个功能：一是在地图上标出每个湖泊的中心位置；二是在用户点击某个湖泊标记时，弹出其对应的表面积变化趋势图。实现这一功能的核心代码如下：

import folium
import plotly.express as px
import os

# 定义函数：绘制时间序列并保存为 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}')
    # 设置 X 轴和 Y 轴标签
    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

# 创建地图实例
m = folium.Map(location=[35.5, -119.5], zoom_start=7)

# 遍历湖泊数据，添加带图表弹窗的标记
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')

# 清理生成的临时 HTML 文件
for spot_id in lakes_df['ID']:
    os.remove(f'tmp_{spot_id}.html')

github.com/stac-utils/pystac-client/blob/main/docs/quickstart.rst

上述脚本中，首先定义了一个绘图函数，用于筛选特定湖泊的数据，使用 Plotly 绘制折线图，并将其导出为独立的 HTML 文件。随后，通过 Folium 初始化一个交互式地图。程序遍历湖泊信息表（lakes_df），为每一个湖泊的地理中心点创建标记，并将之前生成的图表嵌入弹窗中。最终，完整的地图被保存为名为 'map_with_timeseries.html' 的文件。任务完成后，系统会自动删除所有中间生成的临时 HTML 图表文件，以保持环境整洁。

完成以上步骤后，打开项目目录下生成的 HTML 文件，即可看到各湖泊在地图上的分布情况，每个标记代表一个感兴趣区域（AOI）——在此案例中即为湖泊的中心坐标。

www.element84.com/earth-search/examples/

当点击任意一个湖泊标记时，便会弹出该湖泊在 2024 年内由 Sentinel-2 影像计算得出的水面面积变化趋势图。例如，下图展示了编号为 #5 的湖泊的时间序列动态：

2024年湖#7表面面积的时间序列数据由Sentinel-2影像计算得出，图像由作者提供。

towardsdatascience.com/create-an-interactive-map-to-display-time-series-of-satellite-imagery-e9346e165e27

在遥感与地理空间数据分析领域，新技术和工具正以惊人的速度涌现。几乎每个月都有新的Python包或库被开发出来，旨在更高效地提取、处理、可视化大规模数据。尽管工具有所进步，但实际应用中仍面临两大核心挑战：其一是分析人员需具备足够的专业知识，以确保从海量（GB乃至TB级别）遥感数据中提取的信息准确可靠；其二是如何构建一个合理的系统架构，将多个功能库有机整合，从而生成真正有意义且可投入实际使用的成果。

本文通过具体案例强调了图像处理过程中看似微小的操作失误，可能引发最终结果中的显著偏差。同时，在可视化环节，展示了如何将Folium、Plotly等交互式地图工具与最新的图像裁剪API相结合，搭建出一套可用于持续监测地表变化的实用框架。这种集成方式为环境监测、灾害响应等多种应用场景提供了强有力的技术支持。

此外，以下资源链接提供了相关主题的深入探讨：

• 从太空观测风暴：利用Python脚本生成震撼视觉效果
• 一段长Python代码教你如何用地理空间大数据制作短视频
• 卫星如何探测隐匿的熔岩流与活跃野火？——基于Python的技术解析
• 在Google Colab中结合Python与Folium库创建交互式地图的完整指南

github.com/stac-utils/pystac-client/blob/main/docs/quickstart.rst

www.element84.com/earth-search/examples/

参考文献
Copernicus Sentinel 数据 [2024] — 提供Sentinel系列卫星数据
Copernicus 服务信息 [2024] — 关于Copernicus各项服务的官方说明

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：时间序列 交互式 Interactive time-series collections