简介面向地理信息、遥感与生态环境等领域研究者这里提供的是2015年中国区域1km分辨率年度植被指数NDVI空间分布数据集。源数据取自NASA MODIS MOD13A3月合成产品经提取子数据集、逐月拼接、投影栅格、单位换算与全国裁剪并利用最大合成法生成年度NDVI可有效降低云、气溶胶及季节波动干扰。数据采用Albers等积投影WGS84空间分辨率1km时间分辨率为年适合用于植被覆盖度估算、区域生态质量评价、土地利用变化及长时间序列遥感对比等场景。资源包共5个文件压缩后约20.01MB核心是GeoTIFF栅格文件另附TFW定位文件、XML元数据及TXT数据说明便于核对投影参数、数据来源及预处理流程。目前已有312人学习下载能为需要现成全国尺度NDVI底图、希望跳过复杂预处理的研究者、学生或规划人员节省大量制图与分析准备时间。1. 2015年中国 1km NDVI 空间分布数据集先看清你拿到的是一份什么底图做生态评估、植被覆盖度反演或国土空间分析的人大概率绕不开一个需求要一份覆盖全国、时间对齐、分辨率够用又不至于把磁盘塞满的 NDVI 数据。MODIS 2015 年中国 1km 植被指数NDVI空间分布数据集解决的就是这个诉求。它把 MODIS 的 MOD13A2 产品按中国行政边界做了裁剪和重投影形成逐年可比的 NDVI 基底分辨率 1km适合全国及省级尺度的趋势分析、物候提取和 FVC 反演。一个反直觉的事实是1km 看着粗糙但在全国尺度上它比 250m 的 MOD13Q1 少一半以上的存储和运算量又比 10km 的 GIMMS 更能看清地形和土地利用边界。新手可以直接把它当底图用熟手需要关心的是数据背后的合成口径、质量控制层和投影细节。2. 数据源与合成逻辑MOD13A2 的 16 天合成、1km 网格和 QA 层2.1 为什么是 MOD13A2三款 MODIS NDVI 产品的选型对比MODIS 植被指数产品家族里和 NDVI 直接相关的主要有三款MOD13Q1250m16 天合成、MOD13A1500m16 天合成、MOD13A21km16 天合成。标题里写的是“1km 植被指数NDVI空间分布数据集”对应到 MODIS 标准产品就是 MOD13A2。选它而不是另外两款核心原因是研究尺度决定的。做全国尺度的植被动态监测250m 数据单景文件动辄几百 MB全国拼接后处理效率很低10km 的 GIMMS 或 AVHRR 数据虽然轻但空间细节丢失严重连秦岭-淮河这样的重要地理分界线都模糊掉了。1km 是“够看清省级差异、又跑得动”的折中选择。另一个实际考量是时间跨度MOD13A2 从 2000 年持续至今2015 年处于数据质量最稳定的阶段Terra 卫星的平台衰减问题在 2015 年后的定量研究中处理方案已经成熟。如果你做的是长时序趋势分析2015 年这个时间点刚好可以和前后十年无缝衔接。选型时还要注意MOD13A2 的 NDVI 是大气校正后的地表反射率计算而来不是 TOA 反射率直接算的。这意味着它已经考虑了臭氧、气溶胶、水汽的干扰你拿到手后不应该再做一次大气校正否则等于二次处理反而引入误差。部分研究者会在此基础上做 Savitzky-Golay 滤波平滑那是后处理不是预处理。2.2 16 天合成的时间口径2015 年到底对应哪几个合成期“2015 年 NDVI 数据集”这句话听起来简单但时间口径直接决定你后续能不能做物候分析。MODIS 16 天合成产品的年周期是一年 23 个合成期第 1 天到第 353 天附近2015 年对应的就是 2015 年第 1 期到 2015 年第 23 期。常见做法是按“年内合成期序号”组织数据比如 2015_001、2015_017、2015_033……每期之间间隔 16 天但每个合成期的窗口并不是固定从第 1 天开始的。以 2015 年为例第 1 期覆盖 1 月 1 日到 1 月 16 日第 2 期覆盖 1 月 17 日到 2 月 1 日以此类推。合成算法是“最大值合成”也就是在这 16 天窗口内每个像元取 NDVI 最大的那次观测。这个算法的设计初衷是云、雪、阴影都会让 NDVI 偏低取最大值能优先保留晴空观测。这里有一个容易被忽略的坑最大值合成会系统性偏高在植被生长季表现不明显但在荒漠、裸土和冬季积雪区噪声会被放大。所以你在做年度统计时如果直接拿 23 期的最大值作为“年最大 NDVI”结果会和生长季峰值 NDVI 有偏差。更稳的做法是先做 QA 过滤再用生长季5 月到 9 月内的最大值。2.3 1km 背后的投影和网格正弦投影为什么容易让人翻车MOD13A2 原始产品的投影是正弦投影Sinusoidal全球按 10 度经度分带切分为瓦片tile每个瓦片的网格是标准的 1km 近似分辨率。这个投影在赤道附近变形小但到了中国北方像元实际地面尺寸会在东西方向被拉伸直接拿去算面积或和矢量边界叠加会出现明显的错位。我在实际处理中见过不少翻车案例有人直接用 HDF 原始投影做中国区域裁剪出来的边界歪歪扭扭和行政区矢量叠不上还以为是数据有问题。问题不在数据在于投影没统一。做中国区域分析两种常见选择是一是转成 WGS84 经纬度坐标分辨率约 0.0083333333 度约 1km适合和全球尺度数据对比二是转成 Albers 等积投影中央经线 105°E标准纬线 25°N 和 47°N适合做面积统计和省级制图。转投影必然伴随重采样这会带来信息损失。1km 数据重采样到 0.0083333333 度像元大小几乎不变损失很小但如果你转成 Albers 投影且目标分辨率设置不当可能会出现像元拉伸。常见做法是把 Albers 分辨率设置为 1000m而不是 0.01 度之类的混合单位避免网格错位。3. 从 HDF 到全国 NDVI 栅格拼接、重投影、裁剪与 QA 过滤3.1 先摸清 HDF 里的数据层用 gdalinfo 看子数据集结构MOD13A2 的原始文件是 HDF4 格式一个文件里包含多个科学数据集SDS包括 1km NDVI、EVI、像元可靠性Pixel Reliability、质量详情、合成天Composite Day of Year等。用 GDAL 读取时你不能直接把 .hdf 当作一个普通栅格处理必须先列出子数据集。gdalinfo MOD13A2.A2015001.h23v04.006.2015312234157.hdf运行后会在输出中看到类似这样的子数据集路径SUBDATASET_1_NAMEHDF4_EOS:EOS_GRID:MOD13A2.A2015001.h23v04.006.2015312234157.hdf:MODIS_Grid_16DAY_1km_VI:NDVI SUBDATASET_2_NAMEHDF4_EOS:EOS_GRID:MOD13A2.A2015001.h23v04.006.2015312234157.hdf:MODIS_Grid_16DAY_1km_VI:EVI SUBDATASET_5_NAMEHDF4_EOS:EOS_GRID:MOD13A2.A2015001.h23v04.006.2015312234157.hdf:MODIS_Grid_16DAY_1km_VI:pixel_reliability这段命令的作用是查看 HDF 的内部结构确认 NDVI 层和 QA 层的子数据集名称。注意 MOD13A2 的 NDVI 是 int16 类型存储的有效值范围是 -2000 到 10000需要乘以 0.0001 缩放因子才能得到实际 NDVI 值范围 -0.2 到 1.0。很多人在这一步没做缩放直接把整数当 NDVI 用结果算出来的 FVC 全部异常。这也是新手最常犯的错误之一。3.2 统一网格的拼接与重投影gdalwarp 关键参数中国区域在 MODIS 正弦投影下覆盖了大约 15 个瓦片tile大致从 h23v03 到 h28v06 的范围。把这些瓦片拼接成一个全国栅格同时完成重投影最稳的工具是 GDAL 自带的 gdalwarp。示例命令如下gdalwarp -t_srs EPSG:4326 -tr 0.0083333333 0.0083333333 -r cubic \ -of GTiff -co COMPRESSLZW -ot Float32 \ HDF4_EOS:EOS_GRID:...h23v04...hdf:MODIS_Grid_16DAY_1km_VI:NDVI \ HDF4_EOS:EOS_GRID:...h24v04...hdf:MODIS_Grid_16DAY_1km_VI:NDVI \ ... (其余瓦片省略) ... \ china_ndvi_2015_001.tif参数说明-t_srs EPSG:4326把正弦投影转成 WGS84 经纬度-tr 0.0083333333 0.0083333333设置输出像元大小约 1km-r cubic使用三次卷积重采样对 NDVI 这种连续变量比 nearest 平滑但不会像 bilinear 那样过度模糊-ot Float32保留浮点精度但此时数据还是缩放前的整数后面需要用gdal_calc.py或 Python 做缩放。注意瓦片文件较多时在命令行写全路径容易出错我一般先在脚本里把文件列表生成好再拼接成 gdalwarp 的输入串。这个命令完成两件事把多个瓦片拼到一张图里同时把投影转成目标坐标系。如果遇到接缝明显多半是-r选成了 nearest或者不同瓦片之间存在数据重叠但没有做无缝镶嵌。GDAL 在处理重叠区域时会默认取第一个文件的值所以最好在拼接前把瓦片按空间顺序排列好否则边界处可能出现锯齿。3.3 QA 过滤与边界裁剪让雪、云和境外瓦片退场拼接重投影得到的 NDVI 栅格还不能直接用因为云、雪、阴影对应的像元虽然在 NDVI 合成时已经被抑制但部分噪声仍然存在。MOD13A2 的 pixel_reliability 层用 0 到 3 四个等级标记像元质量0 表示理想质量1 表示可用2 表示雪/冰3 表示云/阴影。你需要把这层也按同样的方式拼接重投影然后用来过滤 NDVI。# 1. 将 QA 层拼接重投影保持与 NDVI 完全相同的网格 gdalwarp -t_srs EPSG:4326 -tr 0.0083333333 0.0083333333 -r nearest \ -of GTiff -ot Byte \ HDF4_EOS:EOS_GRID:...h23v04...hdf:MODIS_Grid_16DAY_1km_VI:pixel_reliability \ ... \ china_qa_2015_001.tif # 2. 用 gdal_calc 把质量标记为 0 或 1 的像元保留其余设为无效 gdal_calc.py -A china_ndvi_2015_001.tif -B china_qa_2015_001.tif \ --outfile china_ndvi_2015_001_clean.tif \ --calcA*(B1) \ --NoDataValue-9999这段命令的关键点有两个QA 层重采样必须用-r nearest因为质量标记是类别变量用 cubic 会插值出 0.5、1.5 这种不存在的质量等级而 NDVI 层可以用 cubic因为它本身是连续变量。gdal_calc.py里的B1生成一个布尔掩膜乘在 NDVI 上质量差的位置变成 0但 0 和云区 NDVI 值混在一起会污染统计所以后面要再设置 NoDataValue-9999并在统计时忽略它。裁剪行政边界这一步常见做法是用标准的中国国界矢量数据转成和 NDVI 栅格一致的投影后用 gdalwarp 的-cutline参数完成。注意-cutline要求矢量属性里有一个裁剪属性字段默认用-crop_to_cutline让输出范围贴合边界。如果边界外出现黑色条带多半是 NoDataValue 没有统一检查两个文件的 NoData 设置即可。4. 从 NDVI 计算 FVC像元二分模型的参数标定与批量实现4.1 像元二分模型公式与端元取值的两种路线NDVI 本身是一个无量纲的绿色度指标但很多实际项目需要的是植被覆盖度FVC。这也是“NDVI 计算 FVC”成为高频检索词的原因。最常见的转换方法是像元二分模型FVC (NDVI - NDVI_soil) / (NDVI_veg - NDVI_soil)其中 NDVI_soil 是纯裸土像元的 NDVI 值NDVI_veg 是纯植被像元的 NDVI 值。理论上这两个端元需要实测标定但实际操作中几乎没有人会去逐像元测量常见做法是取研究区内 NDVI 累积分布的低分位数如 5%和高分位数如 95%作为端元。更精细的做法是分土地利用类型分别标定耕地、林地、草地各自有不同的 NDVI_soil 和 NDVI_veg混在一起用全局分位数会把农田的 FVC 拉低、把荒漠的 FVC 拉高。端元取值的差异对 FVC 结果影响极大。我见过同一景影像有人取 1% 和 99% 分位数有人取 10% 和 90%两者算出来的平均 FVC 能差出 0.15 以上。这不是模型本身的错而是端元定义没有统一。在做区域对比时务必固定同一套分位参数否则结论不可比。4.2 批量计算 FVC 的 Python 流程如果你的数据已经完成了第 3 章的预处理得到的是逐年各期 NDVI 栅格接下来就可以批量计算 FVC。下面是一个用 Python 处理全年 23 期数据的流程骨架依赖 rasterio 和 numpyimport rasterio import numpy as np def ndvi_to_fvc(ndvi_path, q_soil0.05, q_veg0.95, nodata-9999): 将 NDVI 栅格转为 FVC 栅格端元取分位数 with rasterio.open(ndvi_path) as src: ndvi src.read(1).astype(np.float32) profile src.profile.copy() ndvi[ndvi nodata] np.nan # 用有效像元的分位数标定端元 valid ndvi[~np.isnan(ndvi)] ndvi_soil np.nanpercentile(valid, q_soil * 100) ndvi_veg np.nanpercentile(valid, q_veg * 100) # 像元二分模型 fvc (ndvi - ndvi_soil) / (ndvi_veg - ndvi_soil) fvc np.clip(fvc, 0.0, 1.0) fvc[np.isnan(fvc)] nodata profile.update(dtyperasterio.float32, nodatanodata) out_path ndvi_path.replace(_clean.tif, _fvc.tif) with rasterio.open(out_path, w, **profile) as dst: dst.write(fvc, 1) return out_path, ndvi_soil, ndvi_veg逻辑说明q_soil和q_veg分别控制裸土和纯植被端元的分位数默认 5% 和 95%可视研究区调整。np.nanpercentile计算端元时忽略了 NoData避免无效值拖累分位数。np.clip把 FVC 限制在 0 到 1 之间防止端元取值异常导致负值或大于 1。需要注意的是这个函数每次调用都会重新计算端元如果你希望全年使用统一端元就把端元计算单独提出来先遍历全年数据算好全局分位数再逐期带入函数。这里的“全年统一端元”非常关键。分 23 期单独取分位数冬季 NDVI 整体偏低端元也会偏低算出来的 FVC 冬季偏高生长季反而偏低完全违背物理意义。正确做法是把全年有效像元合并取生长季内的累积分位数作为全年固定端元。4.3 NDVI_soil 与 NDVI_veg 的取值边界与失败信号端元取值不是越极端越好。NDVI_soil 取 1% 分位数时把水体、云残影、阴影都包含了进去这些像元的 NDVI 可能为负值或接近 0导致 FVC 整体偏高。NDVI_veg 取 99% 分位数时可能只有极少数的密林像元满足造成 FVC 饱和大片中高覆盖度区域全部变成 0.9 以上失去区分度。我的建议是先做诊断输出全年 NDVI 直方图找到两个明显的谷底。裸土在干旱区的 NDVI 通常在 0.05 到 0.15 之间密林通常在 0.7 到 0.85 之间。如果计算出的 FVC 在某个区域出现大面积 0 或 1不是模型错了而是端元取到了异常值。失败信号还包括FVC 均值超过 0.85不太可能出现在全国尺度、FVC 标准差小于 0.05说明端元区间过窄、同一像元在相邻两期之间 FVC 跳变超过 0.3可能是残留噪声或端元不一致。另外FVC 计算对 1km 数据有一个天然限制混合像元问题。1km 像元里通常同时包含裸土和植被NDVI 本身是两者的线性加权合成结果所以像元二分模型在 1km 尺度上是近似估计不确定性在高异质区域会放大。如果你做的是县域或更小范围的研究建议换 MOD13Q1 的 250m 数据。5. 避坑手册五个把 NDVI 用翻车的常见问题5.1 冬季大面积低值不是数据错了是雪没滤干净现象1 月、2 月的 NDVI 栅格里东北、新疆出现大片连续的低值区域数值在 0 到 0.1 之间甚至出现负值和夏季影像对比像“数据空洞”。原因最大值合成能滤掉大部分云但雪和冰的表面反射率在近红外波段显著低于植被导致 NDVI 骤降。如果你没有过滤 pixel_reliability 层里标记为 2雪/冰的像元这些雪区就会保留下来。解决在第 3.3 节的 QA 过滤步骤中把条件从B1改成B0或拆成两步先剔除雪冰像元再剔除云层像元。如果你只有 NDVI 栅格没有 QA 层一个退而求其次的办法是把 NDVI 小于 0.05 且对应时间段为冬季的像元直接设为无效但这种方法要谨慎荒漠地区夏季 NDVI 也可能低于 0.05。5.2 瓦片接缝明显重采样方法和参考网格不一致现象拼接后的全国影像上能看到清晰的瓦片边界尤其是山地和农田区域一条直线两侧的 NDVI 差异明显。原因两种可能性一是gdalwarp的-r参数选了 nearest导致拼接边界处像元取值跳跃二是不同瓦片下载时间不一致虽然同为 2015 年但部分瓦片可能混入了其他合成期的数据。解决重采样改-r cubic并在拼接前用脚本检查每个瓦片的合成期是否一致读取 HDF 文件名里的第 3 段日期码如 A2015001 表示 2015 年第 1 期。若已拼接完成才发现接缝可以用gdal_merge.py配合羽化参数重新处理但效率远低于源头上统一。5.3 中国边界内出现空白或境外数据残留现象裁剪后国界线内出现零星空洞或者边界外有大片非 NoData 的植被值通常是邻国领土上的真实 NDVI。原因边界矢量与栅格投影不一致或矢量边界精度比 1km 像元细得多裁剪时落在像元边缘导致取舍差异。对于境外数据残留纯粹是没有用-crop_to_cutline裁剪。解决裁剪前用ogr2ogr把矢量重投影到与 NDVI 栅格完全一致的坐标系裁剪时用gdalwarp -cutline china.shp -crop_to_cutline -dstnodata -9999。检查输出栅格的统计直方图如果边界外还有大量有效值说明-crop_to_cutline参数没有生效确认 GDAL 版本不要太旧2.2 以下版本对该参数的支持有缺陷。5.4 1km 混合像元让城市、水体、农田齐齐失真现象城市区域 NDVI 普遍低于 0.15但城市公园或大型绿地所在的像元却出现跳变水体边缘出现 NDVI 正值异常。原因1km 像元尺度过大城市中建筑、道路、植被混合在一个像元里NDVI 被建筑反射率拉低水体的近红外反射极低但 1km 像元混合了岸边土壤或水生植被导致 NDVI 虚高。解决凡是涉及城市或水体的研究不要用 1km 数据做像元级分析改用 250m 或 30m 数据。如果只能用 1km至少要做掩膜用水体掩膜和城市不透水面掩膜把这些区域剔除再统计 FVC。掩膜的阈值建议基于同期 Landsat 分类结果提取而不是直接用 NDVI 阈值否则会形成循环论证。5.5 FVC 端元用全局分位数半干旱区直接饱和现象计算 FVC 时内蒙古、甘肃等半干旱区大片像元 FVC 等于 0而云南、广西等湿润区大片像元 FVC 等于 1中间梯度丢失。原因全国尺度的 NDVI 范围极大从荒漠的 -0.05 到热带雨林的 0.9。用全国统一分位数取端元半干旱区的 NDVI 本身就接近全国的低分位数导致该区域的 FVC 全被压到 0。解决分生态区或分省计算端元。常见的分法是按植被地带性划分东部季风区、西北干旱区、青藏高寒区各算一套端元。这样每个区域的 NDVI 动态范围都能充分展开FVC 的区分度显著提升。代价是区域之间的 FVC 统计口径不再完全一致所以如果你要做全国排名或对比需要说明这一方法学选择。6. 验证与进阶这套 NDVI 数据到底值不值得信6.1 三个可执行的验证方法拿到或生产出 2015 年全国 1km NDVI 栅格后先别急着进模型做三个低成本验证能避免后续分析白做。第一个是时间剖面验证在栅格上随机抽 20 到 30 个点覆盖林地、农田、草地、荒漠四类地表提取全年 23 期 NDVI 值画出时间序列曲线。正常曲线应该是一个平滑的单峰或双峰形态峰值在 7 月到 8 月如果出现锯齿状抖动或冬季异常凸起说明 QA 过滤不彻底回炉重做。第二个是和实测站点验证从中国气象数据网或生态站获取 2015 年的植被观测记录如地上生物量、覆盖度实测值和对应像元位置的生长季最大 FVC 做线性回归。R 方大于 0.6 基本可用小于 0.4 就要检查端元取值和混合像元问题。注意站点坐标要精确到秒级1km 像元范围内如果有多个站点取平均值更稳。第三个是高分辨率交叉验证用 2015 年夏季一景无云的 Landsat 8 NDVI30m聚合到 1km 尺度和 MODIS NDVI 做逐像元对比。两者的差值在裸土区通常在 ±0.05 以内在异质区域可能到 ±0.15。这个验证能帮你量化 1km 数据在目标区域的不确定性写论文时可以直接引用。6.2 让数据发挥更大价值的两个进阶方向如果你不满足于把 NDVI 当静态底图2015 年这套数据还可以做两件更有价值的事一是物候提取用 23 期 NDVI 曲线拟合植被生长季的起始期和结束期识别中国植被物候的南北梯度二是和 2000 年、2010 年、2020 年的同期数据做叠加计算 15 年 NDVI 趋势斜率识别生态恢复区和退化区。这两个方向的关键都是保持数据口径一致包括端元参数、QA 过滤强度、重采样方法任何一步变了趋势分析的结果就不可信。我自己的习惯是每套数据到手先建一个“数据卡片”记录投影、分辨率、合成期范围、QA 过滤规则、端元取值。这个习惯已经帮我避免了至少三次返工。数据卡片的成本只有几分钟但能让你在三个月后回头时清楚知道每个栅格是怎么生成的。希望这整套流程能帮到你让你在这条路上少踩几个我踩过的坑。本文还有配套的精品资源点击获取