利用遥感和 Python 追踪休斯顿的绿化覆盖率（基于卫星 NDVI）
变化常常以我们难以察觉的方式发生——除非我们退得足够远。在地球上空，卫星观测着城市的生长、变化和呼吸。得益于欧洲航天局的哨兵二号卫星和免费提供的哥白尼数据，我们现在可以逐像素地追踪城市绿化覆盖的细微变化。
这篇文章探讨的是使用归一化植被指数 (NDVI) 观察休斯顿 2021 年植被状况时会发生什么。NDVI 是光合作用的指标，提取自卫星图像中的红光和近红外波段。植物生长时会反射近红外光并吸收红光。NDVI 将这种对比度转换为 -1 到 +1 之间的值。区域越绿，NDVI 越接近 1。
以下是 2021 年休斯顿的月度数据：
由于完全被云层覆盖，因此五月被排除在外。
你所看到的是休斯顿的绿意盎然，以及这片绿意栖息的土地。公园、河流、未开发的土地——它们在浅绿和浅绿色的色调中格外醒目。树木更茂密的街区则显得不那么红。正如你所预料的那样，高速公路和工业区几乎不会随着季节变化而改变。
这里最震撼的不仅仅是色彩，还有观察的一致性。这张 GIF 由 12 幅 Sentinel-2 图像组成，2021 年每个月一张。每幅图像都是休斯顿核心区的无云快照，使用 NDVI 方程渲染：
NDVI = (NIR - Red) / (NIR + Red)
具体来说，就是用 Sentinel-2 2A 级影像，将第 8 波段减去第 4 波段，再除以它们的和。我们通过 openEO 提取了这些数据，这是一个 API 层，它允许您访问哥白尼数据空间生态系统 (CDSE) 基础设施，而无需在本地存储原始影像。
每一帧都蕴含着某种意义。一月和二月，休斯顿的植被处于休眠状态——一片片红色斑块占据主导地位。但到了三月，万物开始绽放。这座城市焕发出新的生机。我们看到公园里生机勃勃。夏季的绿意一直持续到九月下旬。然后，逐渐褪色的景象又回来了。
NDVI 为何如此有用？
NDVI 是一个定量信号。研究人员通常用它来测量作物健康状况并预测产量、监测森林砍伐和土地利用变化，以及监测城市绿化计划。
在休斯顿这样的地方，开发向外延伸到低密度绿地区域，它是一种追踪我们经常忽略的东西的工具：我们建造的东西。
Sentinel-2 和哥白尼
获取遥感数据的方法有很多。Sentinel-2 非常棒，因为它可以免费获得 10 米分辨率的数据，并且每 5 天提供一次覆盖全球的数据。它提供不同的波段（红光、近红外、短波红外），但比私人服务提供的波段要少。
Copernicus openEO 让您无需下载大文件即可处理影像。因此，这类遥感技术的入门门槛很低。您只需要 Python、一个想法和一些耐心。
接下来会发生什么
这个动画只是一个开始。我正在制作其他示例，展示遥感技术可以实现的更多酷炫功能。
from pathlib import Pathimport openeoimport xarray as xrimport matplotlib.pyplot as pltimport imageioimport numpy as npimport osfrom datetime import datetime, timedelta# Setup foldersbase_path = Path(\"results\")base_path.mkdir(exist_ok=True)frames_dir = base_path / \"frames\"frames_dir.mkdir(exist_ok=True)# Connect and authenticateeoconn = openeo.connect(\"openeo.dataspace.copernicus.eu\")eoconn.authenticate_oidc()# Define Houston bounding box (smaller for speed)bbox = [-95.40, 29.70, -95.30, 29.80] # Downtown Houston# Define year and monthly date ranges (10-day windows)year = 2021intervals = [ ( datetime(year, m, 1).strftime(\'%Y-%m-%d\'), (datetime(year, m, 10)).strftime(\'%Y-%m-%d\') ) for m in range(1, 13)]# Process and download each NDVI tileall_ndvi_slices = []for i, (start, end) in enumerate(intervals): print(f\"\\n Requesting NDVI for {start} to {end}\") # Create NDVI cube cube = eoconn.load_collection( \"SENTINEL2_L2A\", temporal_extent=[start, end], spatial_extent=dict(zip([\"west\", \"south\", \"east\", \"north\"], bbox)), bands=[\"B04\", \"B08\"] ).ndvi(red=\"B04\", nir=\"B08\") # Submit and track job job = cube.create_job(title=f\"NDVI_{start}\", out_format=\"NetCDF\") job.start_and_wait() # Download result outfile = base_path / f\"ndvi_{i:02d}.nc\" job.download_result(outfile) # Load NDVI slice and store it ds = xr.open_dataset(outfile) ndvi_slice = ds[\"var\"].isel(t=0) all_ndvi_slices.append((start, ndvi_slice)) ds.close()# Generate animation framesframes = []for i, (date, ndvi) in enumerate(all_ndvi_slices): fig, ax = plt.subplots(figsize=(6, 6)) ax.imshow(ndvi, cmap=\"RdYlGn\", vmin=0, vmax=1) ax.set_title(f\"NDVI – {date}\") ax.axis(\'off\') fname = frames_dir / f\"frame_{i:03d}.png\" plt.tight_layout() plt.savefig(fname) plt.close() frames.append(imageio.imread(fname))# Save animationgif_path = base_path / \"houston_ndvi_2021.gif\"imageio.mimsave(gif_path, frames, duration=1.0)print(f\"\\n Saved NDVI animation: {gif_path}\")
