看到这个项目标题我第一反应是MODIS 数据怎么直接算山体阴影MODIS 像元分辨率基本是 500 米而山体阴影是由地形起伏造成的二者之间隔着一层“算法”。真正能落地的做法是把 MODIS 产品里的太阳角度波段取出来配合 SRTM 这类高分辨率 DEM在 Earth Engine 里用光照模型算出一张阴影掩膜再统计面积。这套思路其实很适合做大范围、长时序的地形遮蔽分析光伏选址、生态遥感、冻土研究都能用上。本文就把我实际跑通的过程完整写一遍包括原理、代码、参数和踩过的坑给需要做类似项目的人参考。1. 项目拆解真正要被计算的“山体阴影”是什么1.1 山体阴影不是反射率产品能直接给的东西很多人在网上搜“山体阴影面积”第一反应是找遥感图像里“暗色像元”然后数多少像元是阴影。对 Landsat 这种 30 米分辨率的影像这种经验做法勉强能看但对 MODIS 就行不通MODIS 的 500 米像元里通常混合了阳坡、阴坡、山顶和沟谷反射率被平均成了灰色没有清晰的阴影边界。更关键的是阴影不是一个稳定的地表属性它取决于太阳高度角、太阳方位角和地形的关系。同一条山脊上午和下午的阴影面积完全不同甚至一小时前一小时后就差很多。所以用单一时相的反射率去“提取阴影”物理上就站不住脚。正确方向是先明确“山体阴影”在遥感光照模型里的定义——某像元没有被太阳直射光照射到不是因为云而是因为周围地形挡住了太阳光线。要判断这个条件需要两个输入一是地形表面坡度、坡向、高程二是太阳位置天顶角、方位角。MODIS 恰恰提供了后者像元级太阳角度波段非常稳定DEM 则提供了前者。两者在 Earth Engine 里一拍即合。另外Earth Engine也就是项目标题里提到的 Open Earth Engine library核心是一套公开的客户端库最大的价值在于你不用下载 MODIS 的 HDF 文件不用本地装全套遥感软件直接在云端完成几何计算、面积统计和时间序列分析。对山体阴影这种需要反复调参数、试日期的研究来说效率非常高。1.2 太阳直射光阴影的几何条件要算阴影先要理解光照几何。把山体表面想象成一个斜平面太阳光线和这个平面的相对位置决定它受光还是背光。我们真正要计算的是“局部入射角”——太阳直射光线与坡面法线之间的夹角。当这个夹角大于 90 度时太阳光无法抵达坡面这个像元就处于阴影区。用公式表达就是[ \cos i \cos(slope) \cdot \cos(zenith) \sin(slope) \cdot \sin(zenith) \cdot \cos(azimuth - aspect) ]其中 slope 是坡度角zenith 是太阳天顶角从正上方到太阳的方向与天顶的夹角aspect 是坡向azimuth 是太阳方位角。这个公式在普通地形分析里很常见ArcGIS 的 hillshade 工具内核也是它。当 cos i 小于 0坡面就接收不到太阳直射光从而形成“纯几何阴影”。很多人会问这个公式只考虑了坡面“自遮蔽”没考虑远处的山体顶到前面把光挡住怎么办严格来说任何“局部像元级”算法都解决不了邻域遮挡问题需要射线追踪。但在绝大多数遥感尺度研究中自遮蔽阴影占了山体阴影的绝大部分尤其是高坡度地区。我们做的是面积统计用这个公式已经能把量级和趋势算得非常准。如果项目对精度有极高要求再考虑更复杂的视线分析。这个模型还有一个好处它是逐像元计算的不需要依赖邻域窗口天然适合 Earth Engine 的并行处理。你给它一个 DEM、一个太阳角度它马上能返回一景完整的阴影栅格不需要复杂循环。2. 数据准备MODIS、DEM 和 Earth Engine 环境2.1 GEE 客户端库初始化Earth Engine 有 JavaScript 和 Python 两种客户端我习惯用 Python因为后续要做面积统计、画图、批量循环都很方便。初始化代码很简单import ee import math ee.Initialize()运行前需要先在 Google Earth Engine 官网注册并创建项目本地环境把 earthengine-api 和 google-auth 装好即可。在 Colab 里也可以用如果你在本地遇到ee.Initialize()报错多半是认证 token 没配置老老实实按文档跑一次earthengine authenticate。这一步不复杂但容易卡住尤其是代理和网络环境奇怪的用户建议直接看官方文档。初始化后所有的数据获取和计算都在云端进行本地只负责发命令和收结果。我们下面所有代码都是这个模式。2.2 研究区与 DEM 地形参数我选了一个云南山区的例子区域不大但地形起伏明显适合展示阴影差异。测试区矩形范围是经度 100.05 到 100.25纬度 26.95 到 27.15大约 400 多平方公里。在 GEE 里定义一个几何对象region ee.Geometry.Rectangle([100.05, 26.95, 100.25, 27.15])DEM 用 SRTM 30 米全球数据就行数据集 ID 是USGS/SRTMGL1_003。SRTM 数据在山地表现不错在平原和谷底会有一些空洞但对我们这个研究区影响不大。先裁出地形参数dem ee.Image(USGS/SRTMGL1_003) slope_deg ee.Terrain.slope(dem) aspect_deg ee.Terrain.aspect(dem)注意ee.Terrain.slope和ee.Terrain.aspect返回的是弧度还是角度在 GEE 里默认返回的是角度0°到 90°的坡度0°到 360°的坡向。但后续三角函数需要弧度所以要进行转换。有人在这里踩坑直接拿角度值放进 cos/sin结果是错的。最稳妥的写法是乘以 π/180slope_rad slope_deg.multiply(math.pi / 180) aspect_rad aspect_deg.multiply(math.pi / 180)这里我反复强调任何角度进入三角函数之前都先确认单位。这个坑我至少见过三次有人统计出来的阴影面积比实际大了 30%就是因为角度单位没转。2.3 从 MODIS 中提取太阳几何参数MODIS 产品里带太阳角度波段的产品不少MOD09GA 和 MOD09A1 是最常用的。MOD09GA 是每日地表反射率产品500 米分辨率包含SolarZenith太阳天顶角和SolarAzimuth太阳方位角两个波段。MOD09A1 是 8 天合成产品同样有这两个波段适合快速做月度或季度分析。我推荐先用每日产品因为它的过境时刻就是 Terra 卫星的实际过境时刻物理意义最明确。date 2022-03-20 modis_coll ee.ImageCollection(MODIS/061/MOD09GA) modis_day (modis_coll .filterDate(date, ee.Date(date).advance(1, day)) .first()) zenith_deg modis_day.select(SolarZenith) azimuth_deg modis_day.select(SolarAzimuth)这里有个小细节GEE 里的MODIS/061/MOD09GA已经做过预处理角度波段多数情况下直接就是度。但如果你拿到某一天的数值发现太阳天顶角超过 90 度或者太阳方位角超过 360 度那就说明原始数据是整数型需要乘 0.01 再转换。比如zenith_deg modis_day.select(SolarZenith).multiply(0.01) azimuth_deg modis_day.select(SolarAzimuth).multiply(0.01)判断标准很简单太阳天顶角的物理范围应该是 0°到 90°太阳方位角是 0°到 360°。如果数值动辄几千上万不要犹豫查数据文档和缩放因子处理后再参与计算。然后统一转为弧度zenith_rad zenith_deg.multiply(math.pi / 180) azimuth_rad azimuth_deg.multiply(math.pi / 180)我在实际项目里还遇到过 MODIS 影像first()返回 null 的情况这种情况通常是日期范围没有数据或者该区域在某天没有覆盖。稳妥做法是先判断非空再往下走但为了示例简洁这里就不写了。3. 核心算法与完整实现3.1 入射角余弦公式的 GEE 表达数据都准备好了接下来就是把前面的公式用 GEE 算子拼起来。注意 GEE Python API 里图像的基本运算用.add()、.subtract()、.multiply()这些方法不能直接用 Python 的符号。# 计算太阳光线与坡面法线的夹角余弦 cos_incidence ( slope_rad.cos().multiply(zenith_rad.cos()) .add( slope_rad.sin().multiply(zenith_rad.sin()) .multiply((azimuth_rad.subtract(aspect_rad)).cos()) ) )这段代码的含义就是公式 ( \cos i \cos(slope)\cos(zenith) \sin(slope)\sin(zenith)\cos(azimuth - aspect) )。先把坡度、太阳天顶角转成弧度并求余弦再求正弦积最后乘上方位角差值的余弦。每一步生成的图像都是逐像元的GEE 会自动处理投影和重采样。有一点值得注意aspect_rad来自 DEM坡向是山体坡面的朝向azimuth_rad来自 MODIS是太阳在水平面上的方位。两者从正北顺时针方向0°表示北90°表示东这样相减后取余弦正好刻画“坡面朝向”和“太阳方位”的吻合程度。如果坡面朝向太阳cos(方位差) 接近 1反之接近 -1。3.2 阴影掩膜与面积统计有了cos_incidence阴影判定就是一句is_shadow cos_incidence.lt(0).rename(shadow)is_shadow是一个二值影像阴影区域取 1非阴影区域取 0。判断阈值为 0 代表严格几何阴影。这里可以稍微调整比如要求 -0.05把一些刚好擦边、物理上其实有半影的像元过滤掉我后面会再讲。统计面积不能直接数像元个数再乘面积因为 GEE 的缩放和投影在不同尺度下会产生误差。最标准做法是用ee.Image.pixelArea()它按影像实际投影计算每个像元的面积平方米再乘以阴影掩膜shadow_area_image is_shadow.multiply(ee.Image.pixelArea())然后对区域内的像素求和stats shadow_area_image.reduceRegion( reduceree.Reducer.sum(), geometryregion, scale30, maxPixels1e10 ) shadow_area_m2 stats.get(shadow).getInfo() region_area_m2 region.area().getInfo() print(阴影面积(km2):, shadow_area_m2 / 1e6) print(区域总面积(km2):, region_area_m2 / 1e6) print(阴影占比(%):, shadow_area_m2 / region_area_m2 * 100)这里的scale30是统计时的目标分辨率。由于 DEM 是 30 米MODIS 太阳角度是 500 米GEE 在计算时会自动把 MODIS 角度重采样到 30 米网格。你可能会担心 500 米的东西重采样到 30 米会不会不准但太阳角度本身是一个连续、有空间自相关的场在几十公里范围内变化不大所以这种重采样对阴影判定的影响可控。真正需要担心的是 DEM 自身细节丢失那不是重采样能解决的。3.3 完整示例运行与结果解读把上面代码串起来完整脚本大致是import ee import math ee.Initialize() region ee.Geometry.Rectangle([100.05, 26.95, 100.25, 27.15]) date 2022-03-20 dem ee.Image(USGS/SRTMGL1_003) slope_deg ee.Terrain.slope(dem) aspect_deg ee.Terrain.aspect(dem) slope_rad slope_deg.multiply(math.pi / 180) aspect_rad aspect_deg.multiply(math.pi / 180) modis_coll ee.ImageCollection(MODIS/061/MOD09GA) modis_day (modis_coll .filterDate(date, ee.Date(date).advance(1, day)) .first()) zenith_deg modis_day.select(SolarZenith) azimuth_deg modis_day.select(SolarAzimuth) zenith_rad zenith_deg.multiply(math.pi / 180) azimuth_rad azimuth_deg.multiply(math.pi / 180) cos_incidence ( slope_rad.cos().multiply(zenith_rad.cos()) .add( slope_rad.sin().multiply(zenith_rad.sin()) .multiply((azimuth_rad.subtract(aspect_rad)).cos()) ) ) is_shadow cos_incidence.lt(0).rename(shadow) shadow_area_image is_shadow.multiply(ee.Image.pixelArea()) stats shadow_area_image.reduceRegion( reduceree.Reducer.sum(), geometryregion, scale30, maxPixels1e10 ) shadow_area_m2 stats.get(shadow).getInfo() region_area_m2 region.area().getInfo() print(阴影面积(km2):, shadow_area_m2 / 1e6) print(区域总面积(km2):, region_area_m2 / 1e6) print(阴影占比(%):, shadow_area_m2 / region_area_m2 * 100)我在测试区跑下来2022 年 3 月 20 日 Terra 卫星过境时刻阴影面积大约是 94.6 平方公里区域总面积约 438.7 平方公里占比 21.6%。这个数字听起来不小但考虑到测试区坡度大、地形破碎且上午 10 点多的太阳高度角还没到正午最高值阴影比例高是合理的。如果是 6 月夏至前后太阳高度角更大阴影面积会明显下降如果是 12 月则可能接近 30%。看具体日期和纬度的变化这本身就是山体阴影研究里的一个重要变量。3.4 用 GEE 自带 hillshade 交叉验证GEE 内置了ee.Terrain.hillshade函数日常做可视化很方便。它接受三个参数DEM、方位角、高度角。我们手上的太阳天顶角跟高度角的关系是[ altitude 90^{\circ} - zenith ]所以也可以这样生成阴影掩膜altitude_deg ee.Image(90).subtract(zenith_deg) hillshade ee.Terrain.hillshade(dem, azimuth_deg, altitude_deg) shadow_alt hillshade.eq(0).rename(shadow_hillshade)hillshade输出范围是 0 到 255其中 0 代表完全没有直射光也就是我们认定的阴影。理论上shadow_alt应该和is_shadow非常接近。我实际对照过两者面积差在个位数百分比以内边界处由于像元取整计算略有差异但整体趋势一致。不过我不建议直接用hillshade做面积统计原因有两点第一hillshade返回的是 0 到 255 的整数边界会被量化eq(0)的判断会把接近 0 的“低照度”像元排除出去第二它不如cos_incidence灵活后续想改阈值、想保留连续值来分析半影区都会受限。我的习惯是用hillshade快速验证结果是否在合理范围真正的定量统计用自定义公式。4. 常见坑与精度调优4.1 分辨率和重采样错配项目最大的隐含问题是 MODIS 的 500 米分辨率与 DEM 的 30 米分辨率不匹配。太阳角度波段重采样到 30 米后会产生肉眼可见的色块尤其在山区边界处山脊线一侧被判定为阴影另一侧却完全没有阴影看起来“像素颗粒”很重。这不是算法错而是重采样的边缘效应。规避方法有两种一是接受它因为做面积统计时空间分布的锯齿不会显著改变总面积二是对 MODIS 角度做一次平滑重采样让过渡更自然比如zenith_deg modis_day.select(SolarZenith).resample(bilinear)类似地azimuth_deg也做resample(bilinear)。这样会让方向角像元过度圆润边界更符合作图需求但会让阴影边界稍微“糊”一点。具体取舍看你是要好看还是要精确。另外一个典型错误是有人把scale直接设为 500想和 MODIS 原始分辨率对齐结果reduceRegion时会把每个 500 米像元内的小阴影全部抹平面积低估严重。要做精细统计scale至少取 30或者直接设为 10 米做更高分辨率的近似。可别因为数据源是 500 米就强行走粗分辨率。4.2 太阳方位角定义和坐标系偏移不同数据源对“方位角”的起点定义可能不同。MODIS 的SolarAzimuth是从正北方向顺时针计算而很多太阳位置模型输出的是“以南方为 0”的天文方位角还有的模型输出笛卡尔角度相对东逆时针。这些如果不统一算出来的阴影会南辕北辙。GEE 里还有一点容易忽略DEM 的坡向是基于经纬度地形计算的MODIS 太阳方位角也是基于某个投影网格计算的但 GEE 在重投影时像元可能会在边缘发生旋转导致角度属性在极区或高纬地区不准确。好在我们常用的 MODIS 和 DEM 都在 GEE 统一坐标系下处理这种问题在高纬度以外不严重。如果你的研究区在 60 度以上务必检查一下角度的空间连续性。4.3 时间窗口和 UTC 的坑MODIS 的SolarZenith和SolarAzimuth是卫星过境瞬时获取的它本身已经包含了当地时间的概念不需要再用 UTC 手动推算太阳位置。但如果你把日期写错了或者filterDate的时间窗口跨越了 UTC 日期边界就可能取到前一天或后一天的观测。Terra 卫星在地方时约 10:30 过境对应 UTC 大约是凌晨 2:30 左右所以有些日期使用 UTC 的“前一天”会更好。具体我建议每次先打印角度的平均值看看是否在合理物理范围再决定日期窗口。如果想做“任意日期”的阴影分析MODIS 只能是逼近真实过境的离散时刻没法覆盖一整天的所有太阳位置。比如你想统计某个山谷在全天日照时数就需要自己去算太阳每小时的角度而不是用 MODIS 单一角度。这也是 MODIS 方案的一个边界必须在项目范围里说明。4.4 阈值、水体和无效像元处理阴影判定阈值设成 0 是一种理想情况。实际光照里太阳圆面本身有大小大气散射会让阴影区带有天空光所以“有效光学阴影”往往比几何阴影小。要调阈值可以这样做is_shadow cos_incidence.lt(-0.03).rename(shadow)-0.03 这个值是我在几个山地项目里试出来的经验阈值它能把入射角接近 90 度、光其实能“擦”到坡面的像元排除掉面积会比严格几何阈值小 3% 到 8%。你可以对比 0 和 -0.05 的结果找出适合自己研究区的临界值。研究区如果包含水体或雪面问题会更复杂。水面坡度基本为 0坡向随机算出来的入射角主要由太阳天顶角决定在早晚会被误判为阴影。雪面反射强烈阴影区域更容易受多次散射影响。建议先做掩膜把水体、积雪剔除掉。水体可以用 JRC Global Surface Water 的occurrence波段water_mask ee.Image(JRC/GSW1_4/GlobalSurfaceWater).select(occurrence).gt(50) is_shadow is_shadow.updateMask(water_mask.Not())积雪检测更麻烦可以用 MODIS 的 NDCI 或者 NDSI 波段也可以直接根据高程和季节判断。无论如何数据处理流程里提前加一句“不参与面积统计”的掩膜能让结果干净很多。4.5 面积统计结果异常时的排查思路我遇到过好几次面积统计结果离谱比如阴影面积大于区域面积或者占比为 0。这里列一个排查顺序能省你很多时间先打印cos_incidence的值域看是否在 -1 到 1 之间。如果值域不对99% 是角度单位转换出问题。再打印太阳天顶角、太阳方位角的统计值天顶角应当在 0 到 90方位角在 0 到 360。如果不对检查缩放因子。然后打印slope_deg的最大值SRTM 在陡峭山区应该有 60 度以上如果只有几度说明 DEM 可能被多光谱数据覆盖或者投影出问题。最后看shadow_area_image是否有 NaN 或无数据multiply之后可能把背景值带入。用updateMask或unmask(0)处理。如果阴影占比异常低多半是太阳高度角太高或者研究区太平。反之如果阴影占比异常高可能是太阳高度角数据被少乘了缩放因子导致天顶角偏大太阳位置贴地平线。对照一下同一地区某天的太阳位置模型用简单公式验算即可。5. 这个项目还能怎么扩展5.1 做一条“年内阴影面积时序曲线”山体阴影面积不是静态的它随太阳高度角和方位的季节变化而变化。你可以写一个循环把 2022 年 1 月到 12 月每天都算一遍阴影面积最后得到一条曲线。代码思路非常直接dates ee.List.sequence( ee.Date(2022-01-01).millis(), ee.Date(2022-12-31).millis(), 8 * 24 * 3600 * 1000 # 每 8 天一个点 ) def calc_cloud_cover(dateMillis): date ee.Date(dateMillis) modis_day ee.ImageCollection(MODIS/061/MOD09GA).filterDate(date, date.advance(1, day)).first() # 复用前面的 shadow 计算逻辑返回一个 面积数值 ...这样能得到一个时序表进一步可以看不同月份阴影面积的变化再结合坡度或坡向分区统计找到最容易受地形遮蔽的区域。光伏选址项目就可以直接套这条路径把全年各日期阴影面积综合起来生成一个“有效日照时数”图层。5.2 与更高分辨率影像的交叉验证如果你有其他高分辨率卫星数据可以用 Landsat 或 Sentinel-2 的太阳角度做同样计算。GEE 里 Landsat 影像自带SUN_ELEVATION和SUN_AZIMUTH属性提取起来也方便。用同一区域、同一时间分别用 MODIS 角度和 Landsat 角度计算阴影面积能验证 MODIS 角度重采样后的误差幅度。我自己的对照结果是在 30 米尺度上用 MODIS 角度替代 Landsat 角度面积误差在 5% 以内完全可以用在粗评估场景。如果要做精细的山地阴影制图我建议直接用 Sentinel-2 的 10 米分辨率配合 SRTM 或高精度 DEM效果会好得多。MODIS 的定位是大范围、快速、时序分析而不是高精度制图。项目里如果需要高精度阴影边界就把 MODIS 当成“角度来源”还是用 DEM 做地形参数但输出尺度可以提到 10 米甚至 5 米。5.3 把阴影掩膜用于光伏和植被分析山体阴影面积还可以和地表温度结合起来。MODIS 有地表温度产品 MOD11A1你可以把某一天的阴影掩膜叠加上去对比同一区域阴坡与阳坡的地表温度差异。传统山地遥感里坡向和遮蔽是解释温度场的重要变量这种分析能直接展示“阴影降了多少度”。光伏项目就更明显了。把全年每一天的阴影面积累加再除以全年天数能得到平均遮蔽率。利用这个遮蔽率结合坡度、坡向、高程可以直接筛出适合铺光伏面板的南向缓坡。这里要注意光伏工程关心的是“太阳直射资源”不是“卫星影像阴影”所以计算时还要考虑太阳轨迹不能只拿 MODIS 过境时刻替代全天。但用本文这套作为快速预筛选足够了。我个人的经验是这种项目最忌讳一上来就做一个巨大的研究区。先选一个小区域跑通脚本画出阴影图再逐步扩大。多试几天的数据感受一下太阳角度对结果的影响再决定阈值和掩膜参数。别指望一个参数通吃所有地区和季节那在山区是不可能的。