简介面向环境遥感、生态评估与地理信息分析人员提供2015年中国区域1km分辨率NDVI栅格数据。原始数据源自NASA MOD13A3月合成产品经提取子数据集、拼接、投影栅格、单位换算、边界裁剪等步骤并采用最大合成法生成年度植被指数可有效反映全年植被生长状况。投影采用Albers等积圆锥WGS84椭球中央经线105°双标准纬线25°/47°空间分辨率1km适合全国尺度植被动态监测、覆盖度估算和生态变化分析。压缩包内共5个文件整体约20MB核心为TIFF格式栅格附有TFW坐标配准文件、XML元数据及数据介绍文本可直接在ArcGIS、QGIS等GIS软件中加载使用也便于批量处理与制图。目前已有311人学习下载适合高校科研团队、生态环境行业用户进行植被遥感分析、课程实习或区域案例研究能够节省数据预处理时间快速获得可用的年度NDVI图层。1. 为什么搞植被监测的人都绕不开这份2015年的NDVI数据集做生态评估也好搞农业估产也罢只要涉及2015年前后中国地表植被状态你几乎一定会撞上这个数据集。它来自MODIS传感器空间分辨率1km覆盖中国全境用NDVI把植被绿度变成一组可比较的数字。但很多人下载后第一眼就懵了文件名是MOD13A2或MOD13A3开头打开是HDF格式里面没有直接可用的栅格必须先处理投影、尺度因子和填充值。这份数据不是不能直接用而是需要一套固定的预处理流程。本文按我实际跑通过的方案把从下载到出图的每一步、参数怎么设、哪些地方最容易翻车全部摊开讲清楚。2. 认识数据先把MODIS NDVI的“脾气”摸清楚2.1 产品编号与时间分辨率MOD13A2还是MOD13A3标题点名了“1km”和“植被指数NDVI”MODIS陆地标准产品里有好几个候选。常见的是MOD13A2这是16天合成的1km NDVI产品全年23期MOD13A3则是月合成产品分辨率同样是1km全年12期。还有MOD13A1那是500m的16天产品各有各的用场但不属于这个1km范畴。我一般先问自己一个需求问题我要的是物候曲线还是只要月度绿度如果做作物生长监测16天合成的MOD13A2更能捕捉到关键生育期转折如果只是做全国植被覆盖度估算MOD13A3直接按自然月合成省一次时间聚合。需要注意MOD13A2和MOD13A3都用正弦投影数据分幅存储。中国全境大概涉及h23v03到h30v06这些图幅具体看目标区域。下载前建一个表格把产品名、日期、行列号、文件URL列出来比临时一个个找方便得多。2015年这个年份没什么特殊但如果是做前后对比需要保证年份内所有期数都齐全别漏掉某个16天时段。提示MOD13A2的“16天”不是自然月第一次用很容易跟月度数据搞混。可以理解为每年23个“合成期”第1期从1月1日开始每16天滚动一次。2.2 投影与坐标正弦投影如何转成WGS84或AlbersMODIS陆地产品的原始投影是正弦投影很多人在这一步翻车。直接把HDF文件拖进ArcGIS或QGIS看到的坐标值常常是几十万到数百万的假坐标不是经纬度。要得到真实地理位置需要做投影变换把正弦投影转成WGS84经纬度或者直接转成中国地区常用的Albers等积投影。我习惯分两步走第一步先做投影转换第二步再裁切中国边界。如果反过来先裁切再投影边界框在正弦投影下往往是个斜矩形容易把边境地区切掉。投影转换时的目标坐标系做面积统计用Albers等积圆锥投影中央经线105°E标准纬线25°N、47°N做经纬度匹配用WGS84。没有统一要求但整条分析链路里只能选一种不能一会儿经纬度一会儿Albers否则后面叠加气象站点、行政区划图时坐标乱套。2.3 数值范围、尺度因子与填充值HDF里真正要用的子数据集叫NDVI存储为16位整数有效取值范围大致是-2000到10000实际NDVI值需要乘以0.0001这个尺度因子。也就是栅格值0对应真实NDVI值0栅格值10000对应1.0。场景中常见的裸土值约0.05到0.2茂密森林可达0.8以上。直接拿原始整数去计算等于没有量纲很多“NDVI是负值”的新手问题其实只是忘记乘尺度因子。填充值也要单独拎出来。无效像元一般用-3000标记云和水体在某些产品里也可能有特定编码。处理顺序必须是先识别无效值并将其设为NaN然后乘以尺度因子如果顺序反了真实-0.3的NDVI和填充值-3000乘完变成-3000掩膜就彻底失效。还有配套的pixel reliability子数据集取值0到40表示好像元1表示边缘像元2代表被云污染3代表被冰雪覆盖4代表无法判断。后面做时间序列清洗时这个波段比NDVI本身还重要。3. 从下载到出图最小可复现的中国区NDVI处理流程3.1 下载前的准备分幅、时间列表与数据整理下载2015年全年中国区NDVI最麻烦的不是数据量大而是文件名对不上。MODIS分幅命名规则里有水平行号h和垂直行号v中国区大致覆盖h23v03到h30v06之间。建议先根据研究范围确定需要哪些图幅不要一整年所有图幅全部下载那样数据量翻几倍却用不上。常见获取渠道是NASA Earthdata需要注册账号国内也有一些数据镜像站提供整理好的数据集。我一般会把下载任务做成一个文本清单每一行对应一个文件。用浏览器批量下载容易断推荐用脚本配合下载工具避免下载了一半没校验导致文件损坏。损坏的HDF文件在读取时会报错提前校验文件大小能筛掉一部分坏文件。3.2 用Python读取HDF并提取NDVI子数据集这里用Python调用GDAL读取HDF4注意MODIS陆地产品是HDF4格式不是HDF5。核心思路是打开文件后通过子数据集名称定位NDVI再读取对应的数值数组。from osgeo import gdal import numpy as np hdf_path MOD13A2.A2015001.h25v05.006.hdf ds gdal.Open(hdf_path) subdatasets ds.GetSubDatasets() for name, desc in subdatasets: print(name, desc)这段代码先打印出HDF文件里所有子数据集执行后会看到NDVI、pixel reliability等子数据集的全路径。看清楚名称后再单独读取NDVI子数据集。ndvi_name HDF4_EOS:EOS_GRID:MOD13A2.A2015001.h25v05.006.hdf:MODIS_Grid_16DAY_1km_VI:NDVI ndvi_ds gdal.Open(ndvi_name) ndvi_array ndvi_ds.ReadAsArray().astype(np.float32) # 将无效值设为NaN ndvi_array[ndvi_array -3000] np.nan # 乘以尺度因子 ndvi_array * 0.0001逻辑说明先按子数据集路径打开栅格读取原始数值后立刻标记填充值-3000为NaN再乘尺度因子0.0001。之所以在转换前先做掩膜是因为填充值本身带符号乘完尺度因子后仍会污染数据范围。如果后续要输出GTiff需要另外写一个函数把投影信息和地理变换参数写进输出文件单纯保存数组会导致空间参考丢失。3.3 投影转换、拼接与裁剪中国边界单幅HDF只是某个分幅块需要把涉及的图幅投影到同一坐标系再拼接裁剪。我惯用的做法是先用gdal.Warp对每一幅做投影转换并输出GTiff然后用gdal_merge.py把多幅图拼起来最后再用gdal.Warp按中国国界矢量裁剪。from osgeo import gdal src gdal.Open(ndvi_name) warp_options gdal.WarpOptions( formatGTiff, dstSRSprojaea lat_125 lat_247 lon_0105 datumWGS84, xRes1000, yRes1000, resampleAlgbilinear, outputTypegdal.GDT_Float32 ) gdal.Warp(ndvi_project.tif, src, optionswarp_options)参数说明dstSRS指定Albers等积投影xRes和yRes表示输出像元大小这里设置为1000米和原始1km分辨率保持一致。resampleAlg选bilinear双线性插值NDVI这种连续变量用双线性比最近邻更平滑但如果后续要做面积统计建议改用cubic或平均方法减少边缘锯齿。做完这一幅后把同一天的所有分幅都投影输出然后执行合并。gdal_merge.py -o ndvi_merged.tif ndvi_project1.tif ndvi_project2.tif这条命令把所有投影后的分幅合并成一个文件。合并后坐标系已经统一接下来加载中国边界矢量做裁剪。注意裁剪时要选与栅格相同的投影否则边界位置对不上。3.4 把16天合成数据整理成月或年NDVI当拿到23期16天NDVI后很多应用并不需要逐期分析而是需要月均值、年最大或季节累计。这里需要明确一个原则MODIS的16天合成本身已经做过一次最大值合成MVC你在做月合成时不能简单把所有16天值取平均那样会把云残渣又带回来。import numpy as np # monthly_arrays: 一个列表包含当月相对的2-3期NDVI二维数组 # 先沿时间维取最大值实现对残留云的进一步压制 monthly_max np.nanmax(np.stack(monthly_arrays, axis0), axis0)逻辑说明对当月包含的2到3期合成结果逐像元取最大值原理是如果某期被云污染NDVI会异常偏低取最大值能选到植被信号最强的那个时刻。这也是MODIS官方在月度产品里采用的核心思路。年最大合成也相同只是把23期全部堆到一起。如果要做年均值我建议先逐期去云确认质量标记之后再做平均否则几个坏像元会把均值拉低一截。4. 预处理避坑这些年我踩过的NDVI数据坑4.1 现象NDVI值出现负值和超过1第一次处理时我把NDVI栅格直接拉伸显示结果发现像元值从零下几千到几万都有。原因是对原始整数没做尺度因子转换还混进了填充值。解决方法是严格按“先掩膜后缩放”的顺序并且把数据范围限定在[-0.2, 1.0]之间超过这个区间的在应用层视为异常。后来我习惯在处理后的每个文件上打印最小值最大值一旦看到负数绝对值特别大立刻回头检查掩膜。4.2 现象时间序列曲线突然掉到0.1以下某像元的NDVI月序列在生长季出现一个明显低值像冬小麦返青时节突然掉到0.1这几乎可以断定是云污染没清理干净。16天合成已经降低了云影响但厚云、冰雪残留仍会漏网。解决方式是读取pixel reliability波段把标记为2或3的像元直接置空再对时间序列做滤波。只靠NDVI自身阈值很难彻底解决因为云在高山区域和低植被区的数值可能重叠。4.3 现象投影转换后图像出现斜向拉伸和裂缝最初做投影转换时我直接用默认的最近邻重采样结果在山区出现明显的锯齿和斜状裂缝。原因是对正弦投影的理解不到位重采样方法又太粗糙。解决方式是改双线性或三次卷积重采样同时把目标像元分辨率设为1000m不要用默认的原生分辨率。另外要注意输出边界范围如果栅格范围设置得太小边缘像元会直接丢失。4.4 现象中国国界附近出现白色条纹裁剪后国界线内侧有一圈无值像元看起来像给中国地图描了个白边。原因是MODIS陆地产品的陆地掩膜本身与国界线不完全一致沿海和边境的相邻像元可能被掩膜标为无效投影后拉成白边。解决方法是在投影转换时设置srcNodata并指定有效数据范围但更彻底的办法是不要过度依赖数据自带掩膜而是用自己准备的边界矢量做裁剪对内部无值像元用邻近像元插值补充。4.5 现象与Landsat或气象站点数据对比时差异大拿1km NDVI与30m Landsat NDVI对比误差往往在0.2以上。原因是空间分辨率不同混合像元效应导致MODIS更平滑。这不是数据错误不能简单地做像素级对比。解决方法是把Landsat NDVI先聚合到1km或者只比较同质性较好的像元比如连片农田内部而不是农田和村庄混在边缘区。只要有分辨率差异就永远不要要求数值完全一致。5. 让数据出活NDVI在农业与生态场景的典型用法5.1 不同作物NDVI曲线怎么用阈值区分作物类型拿到2015年全年NDVI时间序列后很多人的第一反应是计算年累计值但这样做会丢失物候信息。不同作物的NDVI曲线差别非常明显比如冬小麦是典型的双峰或单峰前移峰值出现在4到5月夏玉米峰值出现在7到8月水稻的峰形更宽插秧期有一个明显低谷。要做作物分类重点提取三个指标峰值NDVI、峰值出现时间、生长季NDVI积分。# 假设ndvi_ts是长度为12的月NDVI数组月份索引0-11 peak_value np.nanmax(ndvi_ts) peak_month np.nanargmax(ndvi_ts) 1 # 生长季积分4月到10月的NDVI之和 season_sum np.nansum(ndvi_ts[3:10])参数说明peak_value反映作物最高覆盖度peak_month是峰值出现的月份冬小麦的峰值月一般是4或5夏玉米是7或8。season_sum用于区分一年一熟和一年两熟两熟地区往往有两个峰或更长的生长季。我实际用这三组特征做随机森林分类在黄淮海平原对冬小麦的总体精度在85%左右比直接用最大NDVI好很多。5.2 最大值合成与季节累计参数选16天还是月合成做季节累计时有人直接用MOD13A3月产品累加有人在16天产品上手动聚合。两者结果差别不大但逻辑不同。月产品已经做过MVC做季节累计时直接对月NDVI求和会平滑掉很多异常低值适合大范围生态监测。16天产品做累计前必须完成云污染清除否则某期坏数据会直接进累计值。我一般建议如果目标是全国尺度的植被生产力估算用MOD13A3月产品如果目标是作物单产趋势或物候分析用MOD13A2的16天产品。千万别把两类产品混在同一时间序列里因为它们的合成周期不同插值后会出现系统性偏差。5.3 结合地表温度算植被状态指数地表温度数据可以直接用吗很多做干旱监测的人把“MODIS下载地表温度数据可以直接用吗”挂在嘴边原因是LST产品同样需要尺度因子、投影转换和云掩膜和NDVI的预处理几乎一个套路。植被状态指数通常需要NDVI和LST联合计算比如VCI基于NDVI的多年最小最大归一化TCI基于LST的多年最大最小归一化。vci (ndvi_current - ndvi_min) / (ndvi_max - ndvi_min) tci (lst_max - lst_current) / (lst_max - lst_min)参数说明两个指数都依赖多年基准期所以不能只靠2015年一年数据至少需要2000到2020年一套时间序列来构建最小最大基准。2015年作为待监测年份计算VCI后再与同期降水数据做相关性验证能明显看出农业干旱的空间分布。这里特意强调一点LST数据下载后不能直接读取白天地表温度产品受云影响极大必须配合质量波段使用否则VCI计算结果会被异常高温像元毁掉。6. 最后一个技巧交付前把时间序列“洗”干净在日常项目中我会要求最后交付的NDVI结果通过一道“清洗”流程。注意这里不是简单改数据而是验证后写入元信息让下游同事不会再问奇怪的问题。第一步用pixel reliability波段做硬掩膜把云、雪、坏像元全部设为NaN然后对时序做Savitzky-Golay滤波。SG滤波适合NDVI这种季节性曲线窗口宽度通常取5到7期对于16天产品我常用窗口5、多项式阶数2过宽的窗口会削平峰值过窄又滤不掉噪声。第二步对滤波后的结果重新计算峰值和生长期累计值并把生成的掩膜文件一并交付。掩膜文件可以是一个0/1栅格1表示该像元在整个2015年至少有一期有效观测。训练习惯是这样的每处理一批数据就把最小值、最大值、有效像元占比打印出来存档下次接手同一套流程时能快速发现数据异常。不要只把NDVI扔给同事而是在数据旁边写一个一行描述有效范围、投影、尺度因子、无效值情况。这些细节不值钱但能在团队协作里少几次返工。希望帮到你。本文还有配套的精品资源点击获取
