4. 云阴影几何投影掩膜生成

几何投影的基本关系

云阴影的位置由“云体遮挡直射光”这一几何事实决定：若把云层近似为高度 H 的水平薄层，太阳方位角为 Az（自正北顺时针）、太阳高度角为 El，则云在像元 (x, y) 处的阴影会沿太阳入射方向的反方向落到地面，水平位移量为：

L = H / tan(El)，位移方向为太阳方位角的反方向（即背向太阳）。

在影像坐标系中（x 向东、y 向北），阴影相对云的偏移为：

dx = -L · sin(Az)，dy = -L · cos(Az)

注意符号约定：方位角自正北顺时针度量，因此东向分量为 sin(Az)，北向分量为 cos(Az)；若影像行号向下增长（常见栅格约定），需将 dy 取反后再换算为行偏移。本节统一采用“x 向东、y 向北”的数学坐标，最后再转换到行列号。

从云掩膜到阴影候选区

实际操作中不需要逐像元求交，而是把云掩膜整体平移 (dx, dy) 得到阴影候选区。关键点有三：

• 平移量以像元为单位，需用像元大小 res（米/像元）把 L 换算成像素数：L_px = L / res。
• 平移后超出影像边界的部分应裁剪，边缘处阴影可能落在影像外，属于正常缺失而非错误。
• 云掩膜本身可能包含薄云、碎云，投影前建议先做小半径闭运算，避免阴影区出现大量孔洞。

云高误差对阴影偏移的敏感性

由 L = H / tan(El) 可知，阴影位移对云高是线性敏感的：云高误差 ΔH 引起的位移误差为 ΔL = ΔH / tan(El)。太阳高度角越低，同样的云高误差被放大得越厉害。例如 El = 30° 时 1/tan(El) ≈ 1.73，而 El = 60° 时仅约 0.58。这意味着低太阳高度角条件下，云高估计不准会直接导致阴影掩膜整体错位，与真实阴影的重合度显著下降。

可执行的投影实现

下面给出一个不依赖 GIS 库的纯 Python 实现，输入为云掩膜（二维布尔数组）、太阳方位角、高度角、云高与像元大小，输出阴影候选掩膜。为便于验证，同时给出一个合成小例子。

import numpy as np

def cloud_shadow_projection(cloud_mask, sun_az_deg, sun_el_deg,
                            cloud_height_m, pixel_size_m,
                            row_increases_south=True):
    """
    cloud_mask: 2D bool, True=云
    sun_az_deg: 太阳方位角，自正北顺时针
    sun_el_deg: 太阳高度角
    cloud_height_m: 云高估计值（米）
    pixel_size_m: 像元大小（米/像元）
    row_increases_south: 影像行号是否向南增长（常见栅格为 True）
    返回：阴影候选掩膜（bool），与输入同尺寸
    """
    if sun_el_deg <= 0:
        raise ValueError("太阳高度角必须为正")

    # 水平位移（米）与像素数
    L_m = cloud_height_m / np.tan(np.deg2rad(sun_el_deg))
    L_px = L_m / pixel_size_m

    az = np.deg2rad(sun_az_deg)
    # 数学坐标：x 向东，y 向北；阴影沿太阳反方向
    dx = -L_px * np.sin(az)          # 列方向（东为正）
    dy_north = -L_px * np.cos(az)    # 北向为正

    # 转换为行列偏移
    col_shift = int(round(dx))
    if row_increases_south:
        row_shift = int(round(-dy_north))   # 北向为正 → 行号向南增长需取反
    else:
        row_shift = int(round(dy_north))

    # 平移云掩膜
    shadow = np.zeros_like(cloud_mask, dtype=bool)
    rows, cols = cloud_mask.shape
    r_src0 = max(0, -row_shift); r_src1 = min(rows, rows - row_shift)
    c_src0 = max(0, -col_shift); c_src1 = min(cols, cols - col_shift)
    if r_src1 > r_src0 and c_src1 > c_src0:
        shadow[r_src0 + row_shift:r_src1 + row_shift,
               c_src0 + col_shift:c_src1 + col_shift] = \
            cloud_mask[r_src0:r_src1, c_src0:c_src1]
    return shadow, (row_shift, col_shift), L_px

# 合成示例：一块 5x5 云位于影像中部
cloud = np.zeros((40, 40), dtype=bool)
cloud[15:20, 15:20] = True

shadow, shift, L_px = cloud_shadow_projection(
    cloud, sun_az_deg=135, sun_el_deg=35,
    cloud_height_m=2000, pixel_size_m=10,
    row_increases_south=True)

print("位移像素数 L_px =", round(L_px, 2))
print("行列偏移 (row, col) =", shift)
print("阴影像元数 =", shadow.sum())

示例参数下 L = 2000 / tan(35°) ≈ 2856 m，L_px ≈ 285.6 像素，远超 40×40 的合成影像尺寸，因此平移后阴影会完全移出边界、像元数为 0。这恰好说明：当云高与像元大小不匹配时，必须检查位移是否超出影像范围。实际工程中应把云高、像元大小与影像幅宽一起核对，或改用更小的云高/更大的像元做演示。

云高误差的量化评估

为评估云高误差影响，可对一组候选云高做敏感性扫描，统计阴影掩膜相对基准的位移偏差与重叠度（IoU）：

def iou(a, b):
    inter = np.logical_and(a, b).sum()
    union = np.logical_or(a, b).sum()
    return inter / union if union else 1.0

base_H = 2000
base_shadow, _, _ = cloud_shadow_projection(
    cloud, 135, 35, base_H, 10, True)

for dH in [-500, -250, 0, 250, 500]:
    s, sh, _ = cloud_shadow_projection(
        cloud, 135, 35, base_H + dH, 10, True)
    print(f"ΔH={dH:+5d} m  偏移={sh}  IoU={iou(base_shadow, s):.3f}")

预期现象：ΔH 越大，行列偏移变化越明显，IoU 越低；在低太阳高度角（如 20°）下，同样的 ΔH 会造成更大的偏移差。判断依据是：若业务允许的阴影定位误差为若干像元，则可由 ΔL = ΔH / tan(El) / res 反推可接受的云高误差上限。