Python 追踪休斯顿的绿化覆盖率

这里最震撼的不仅仅是色彩，还有观察的一致性。这张 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 Path
import openeo
import xarray as xr
import matplotlib.pyplot as plt
import imageio
import numpy as np
import os
from datetime import datetime, timedelta

# Setup folders
base_path = Path(“results”)
base_path.mkdir(exist_ok=True)
frames_dir = base_path / “frames”
frames_dir.mkdir(exist_ok=True)

# Connect and authenticate
eoconn = 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 = 2021
intervals = [
(
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 tile
all_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 frames
frames = []
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 animation
gif_path = base_path / “houston_ndvi_2021.gif”
imageio.mimsave(gif_path, frames, duration=1.0)
print(f”\n Saved NDVI animation: {gif_path}”)