简介：提供2010年中国地区1km分辨率NDVI年度空间分布数据集，适合遥感、生态、农业、气候等科研人员与GIS分析用户直接使用，也可作为高校相关课程实习数据。该产品基于NASA MODIS MOD13A3月度植被指数，经过子数据集提取、影像拼接、Albers等积圆锥投影转换、单位换算以及中国区域裁剪，再通过最大值合成法生成全年NDVI栅格，整体已配齐坐标参考与元数据，可无缝嵌入ArcGIS、QGIS等常用平台。压缩包共5个文件，核心为TIFF栅格数据、配套TFW世界文件、两个XML元数据文档及TXT数据介绍，包体约19.96MB，轻量且便于本地存储与共享。目前已有306人学习使用，尤其适合需要现成年度植被指数底图、希望免去复杂遥感预处理流程的学者与开发者。

1. MODIS 2010年中国1km植被指数（NDVI）空间分布数据集：看似现成，实则从分幅到成图每一步都在藏坑

MODIS 2010年中国1km植被指数（NDVI）空间分布数据集，表面上是把分幅产品裁剪后拼起来就能用，可真正动手做的人会发现，光是把几十个HDF文件变成一张全国1km的NDVI栅格，就够耗掉两三个工作日。做植被长时序分析的人最容易在这个年份栽跟头：网上免费能下的中国年NDVI多数是8km的GIMMS，或者是重采样到0.05度的气候网格，跟站点数据一对比尺度就露馅。真正1km的数据必须自己从MODIS处理出来。这份数据集解决的是"2010年中国全境、月度及以上尺度、1km空间粒度"这类需求，省去了分幅下载、投影变换、边界裁剪的重复劳动，但拿到手不等于能直接用：比例因子、无效值、云污染标记都藏在HDF的角落里。以下按我平时做这批数据的处理习惯，把选型、下载、投影裁剪、月度合成和排错串成一条可复现的流程，适合有GIS基础、要做时序统计或生态建模的人。

2. 数据源与产品选型：为什么默认从MOD13A3开始，而不是MOD13Q1

2.1 MODIS NDVI的物理含义与1km尺度由来

MODIS的NDVI不是直接用原始辐亮度算的，而是先对L2反射率做大气校正，再在合成周期里挑选"最干净"的像元。公式是 NDVI = (ρ_nir - ρ_red) / (ρ_nir + ρ_red)，叶绿素对红光吸收强、对近红外反射强，所以健康植被NDVI接近0.8，裸土接近0，水体通常为负。1km这个尺度对应MODIS的MOD13A3产品，比250m的MOD13Q1计算量小，比8km的GIMMS能捕捉到县级边界，做2010年静态分析正合适。

MOD13A3的合成算法不是简单平均，而是最大合成法（MVC）：在合成窗口内保留每个像元NDVI最大值。这套逻辑从AVHRR时代就开始用，目的是消除云和气溶胶的瞬时干扰。MODIS在此基础上增加了逐像元的QA信息，这才是1km数据集真正值钱的地方——只有NDVI值没有质量标记的栅格，拿到手等于一个黑匣子，你不知道哪个像元是被云污染过的。

2.2 产品选型对照：分辨率、合成周期与应用场景的取舍

2010年中国1km植被指数，直接数据源绕不开MOD13A3。和它容易混淆的是MOD13A1、MOD13Q1，以及Aqua卫星的MYD13A3。选错产品会导致后续所有处理链推倒重来，所以先把对照表放出来：

为什么优先MOD13A3？它已经完成了月度合成，下载12幅影像就能做年值。MOD13Q1虽然分辨率更高，但要先从16天合成到月，再做年至，每一步重采样都会损失精度，还要处理250m与1km像元对齐的问题。如果是研究2010年中国全境，我一般只用MOD13A3，Aqua的MYD13A3留作验证。Terra和Aqua的NDVI在相同月份的差异通常在0.05以内，若发现差异过大，优先怀疑云覆盖和分幅拼接顺序，而不是算法本身。

2.3 下载与目录组织：用pyModis批量拉取2010年全年数据

MODIS C6产品的获取主要走NASA LP DAAC，直接在Earthdata Search按时间、分幅筛选也能下，但中国区域至少涉及 h23v05 到 h28v07 的十几个分幅，手工下载一个月点几十个链接，效率太低。常见做法是用 pyModis 这类CMR客户端，它替你处理了Earthdata登录和token。

python modisDownload.py -p MOD13A3 -t h26v05,h26v06,h27v05,h27v06 \ -s 2010-01-01 -e 2010-01-31 -u 你的用户名 -w 你的密码 \ -d /data/modis/2010

逻辑说明：-p指定产品名；-t指定分幅号，中国1km月度产品主要覆盖 h26v05、h26v06、h27v05、h27v06 这几个；-s和-e限定时间窗口；-d是下载目录。pyModis 会先查CMR索引再下载，下载返回的文件名里带实际观测日期，类似MOD13A3.A2010001.h26v05.006.*.hdf，其中2010001表示1月1日，但生产时间戳会变，不能靠硬编码文件名去猜。

下载完成后要检查数量，这一步别跳过：

find /data/modis/2010 -name "MOD13A3.A2010*.hdf" | wc -l

12个月、N个分幅，数量应该是 12 × 分幅数。如果对不上，大概率是某个月的CMR记录缺失或下载被中断，后面所有处理等于在黑洞上盖楼。下载目录最好按产品/年份/分幅划好，因为投影、裁剪和拼接要同时读多个分幅，目录一乱，中间GeoTIFF会连自己也找不到。

3. 从HDF到可用的中国1km栅格：投影、裁剪与单位换算的一体化流程

3.1 HDF4文件结构：NDVI藏在哪个数据层

MOD13A3文件是HDF-EOS格式，GDAL能打开，但要先搞清里面有哪些科学数据集。直接用QGIS拖进去发现是黑白条，就是因为读到了整幅HDF的默认视图而不是NDVI子集。先用pyhdf做一次"解剖"：

from pyhdf.SD import SD, SDC f = SD('MOD13A3.A2010001.h26v05.006.*.hdf', SDC.READ) print(f.datasets().keys()) # 常见输出: ['NDVI', 'EVI', 'DetailedQA', 'pixel_reliability', 'composite_day_of_year', ...] ndvi = f.select('NDVI') raw = ndvi.get() attrs = ndvi.attributes() print(attrs['scale_factor'], attrs['_FillValue'], attrs['valid_range']) # 典型输出: 0.0001 -3000 (-3000, 10000)

逻辑说明：pyhdf 读到的是原始整数，scale_factor=0.0001表示真实NDVI = 原始值 × 0.0001；_FillValue=-3000是无效像元；valid_range给出有效范围。这些元数据在后续单位换算中必须保留，不能等到裁剪后丢掉。这里容易漏：有些教程只把大于0的像元保留，却不管 -3000 和 -2864 这类边缘填充值，最终统计会多出一堆离群值。

3.2 用gdal_translate抽出NDVI子数据集并保持整型

HDF里的投影是MODIS正弦投影，直接拿原始HDF和矢量叠加会对不上。第一步是抽出子数据集并转成GeoTIFF：

gdal_translate -of GTiff -ot Int16 \ HDF4_EOS:EOS_GRID:"MOD13A3.A2010001.h26v05.006.*.hdf":MOD_Grid_monthly_1km_VI:"NDVI" \ /data/work/ndvi_raw_2010_001_h26v05.tif

参数说明：-ot Int16把NDVI保持为整型，避免重投影时插值改变有效值；中间那段HDF4_EOS:EOS_GRID:...是GDAL访问HDF-EOS网格的固定写法，网格名MOD_Grid_monthly_1km_VI在文件metadata里写死。如果报错，先执行gdalinfo MOD13A3.A2010001.h26v05.006.*.hdf查看实际的网格名和子数据集名称。这里用通配符匹配生产时间戳，前提是目录里只有一个该日期文件，否则GDAL会拒绝打开。

3.3 重投影到WGS84并按中国边界裁剪

接下来进行全国范围的投影统一和裁剪。1km数据在中国区域通常用经纬度坐标表达，目标投影直接选EPSG:4326。

gdalwarp -t_srs EPSG:4326 -tr 0.01 0.01 -r near \ -cutline china_boundary.shp -crop_to_cutline -dstnodata -3000 \ /data/work/ndvi_raw_2010_001_h26v05.tif \ /data/work/ndvi_2010_001_h26v05_wgs84.tif

这里至少有三个参数需要每次处理前确认。第一个是-tr 0.01 0.01：0.01度在赤道约1.1km，纬度越高网格越密，对全国范围是常用近似；要严格按1km均匀网格，得用0.008333度，但中国高纬地区文件会变大，实际很少有人这么做。第二个是-r near：默认重采样是cubic，对NDVI这种连续场，cubic容易在云边缘产生振铃，造成 -0.2 的假值，用near能保证不制造新数值。第三个是-crop_to_cutline：如果不指定，输出范围是栅格与cutline面积交集的外接矩形，之后还得用掩膜再裁一次；直接指定后按cutline形状裁剪，代价是边缘像元可能被切掉一点，批量处理时可接受。

如果想要把分幅拼接和裁剪合并，顺序很重要。我的顺序是：先统一投影，再拼接，再裁剪。如果先裁剪再拼接，分幅接缝处会出现重叠区，后续最大合成时选值不统一，接缝痕迹会留在年值图里。

4. 月度合成与年值计算：最大合成法、QA掩膜与堆栈输出

4.1 为什么不能直接用12个月的算数平均做年值

拿到12个月的中国1km NDVI GeoTIFF后，最容易掉进去的坑是直接算年平均。年平均会把冬季低值拉下来，而最大合成保留的是每个月"最像植被"的那天，更适合表达植被覆盖状况。MODIS的MVC算法本身就是在不同观测里挑最大值，我们做年值阶段再次使用MVC，相当于在月度尺度上做了二次选择。这样处理的多云地区年值会偏高，因为每个月只要有一次晴天被选中，全年就是高值。若在研究里声明用的是年最大合成，要在方法里写清楚，否则审稿人会质疑你云污染处理不彻底。

4.2 用Python做年度最大合成：从12个月文件生成年值栅格

以下代码把12个月文件全部读进内存，做逐像元最大值。只需要保证所有文件行列数一致，这在统一投影和裁剪后已经满足。

import numpy as np import glob import rasterio files = sorted(glob.glob('/data/work/ndvi_*_wgs84.tif')) with rasterio.open(files[0]) as src: meta = src.meta.copy() meta.update(dtype='int16', nodata=-3000) # 将12个月叠成三维数组，注意所有文件的行列数必须一致 stack = np.stack([rasterio.open(f).read(1) for f in files], axis=0) valid = stack > -3000 # 只有有效像元参与最大值 any_valid = valid.any(axis=0) annual = np.full((valid.shape[1], valid.shape[2]), -3000, dtype=np.int16) annual[any_valid] = stack[:, any_valid].max(axis=0) with rasterio.open('/data/work/NDVI_2010_annual_max.tif', 'w', **meta) as dst: dst.write(annual, 1)

逻辑说明：stack的形状是 (12, rows, cols)。valid掩膜把所有填充值排除掉，max只在any_valid为True的像元上做，避免全云区被max算成 -3000。annual初始化成 -3000，后续统计时可直接当nodata处理。参数说明：dtype保留 int16，因为还没乘 scale_factor，整数栈更省空间；一旦在年值图上乘了0.0001，就必须转成float32，否则小数部分被截断。

如果内存不够，可以改用rasterio的窗口循环读，分块计算。全国1km大约1500万像元，12层约360MB，绝大多数机器能一次读完。这个环节最容易出现的错误是文件顺序没排对，sorted(glob.glob())按字符串排序，但文件名的生产时间戳并不是观测时间，所以务必先按观测日期重命名文件，或直接把日期提取出来排序。

4.3 QA波段与冰雪污染：pixel reliability的位运算

真正的拦路虎是云标记。MOD13A3的pixel_reliability图层是逐像元的0/1/2/3质量标记，0表示好，1表示可用，2表示云/阴影，3表示冰雪。合成前先用它把2和3剔除，否则冬季冰雪像元会以很高的NDVI破坏年最大值。看一段处理代码：

import rasterio with rasterio.open('/data/work/ndvi_qa_2010_001.tif') as src: qa = src.read(1).astype(np.uint8) reliability = qa & 0b00000011 # 提取最低两位 good = (reliability <= 1)

逻辑说明：pixel reliability的信息被压缩在最低两位，高位还有其他标记。直接对整行QA做位与，只取低两位，就不会把相邻的"雪覆盖"误判成云。注意：如果源文件是MOD13A3的HDF，pixel reliability本身已是一个独立图层，不需要再对DetailedQA做复杂位运算。很多教程推荐"quality=0"严格滤波，会把大量山地数据变成空洞；用<=1能保留必要样本，代价是混入少量云边像元，实际使用时更平衡。

4.4 单位换算与月度堆栈输出

最后把有效NDVI乘以0.0001，输出成float32的月度堆栈，后续时间序列分析就不用再惦记原始编码。这一步我习惯用gdal_calc.py批量处理：

gdal_calc.py -A ndvi_2010_001.tif --outfile=ndvi_f_2010_001.tif \ --calc="A*0.0001" --NoDataValue=-0.3000

参数说明：-0.3000是原始填充值-3000乘以比例因子后的结果。把它设置为nodata后，在Python里用np.nanmask就能直接处理，不用再记忆原始编码。所有12个月都执行一遍后，再做一次年度最大合成，得到的就是真正意义上的浮点NDVI年值图。

5. 常见问题排查：跑2010年数据最容易翻车的5个地方

5.1 整幅影像显示为深黑色，直方图全挤在负值区间

现象：在ArcGIS或QGIS里打开，显示值范围是 -3000 到 1，拉伸直方图还是黑乎乎一片，看不到地表轮廓。原因：没有乘比例因子，或把填充值当成真实值参与统计。NDVI真实值范围应该在 -0.1 到 0.9 之间，若原始值在 -3000 附近成片出现，说明裁剪后无效值没被正确设置成nodata。解决：先建掩膜，再乘scale_factor，最后设置nodata。不要在原文件上直接修改，保留一份Int16原始底图，方便追溯。

5.2 分幅拼接后接缝出现亮线或暗线

现象：两幅相邻分幅的交界处有一条明显折线，在年值图上尤其刺眼。原因：相邻分幅的16天合成时间窗不完全一致，重叠区像元取自不同日期；另外gdal_merge.py默认重叠区用第一个文件的值，简单覆盖而不是融合。解决：先在统一投影网格上做一次"重叠优先"处理，用numpy的maximum在重叠区取两者最大值，而不是依赖gdal_merge的默认行为。代码上就是把分幅先读成两个numpy数组，重叠区做np.maximum(a, b)，非重叠区保持原值。

5.3 沿岸线少了一排像元，岛屿消失

现象：裁剪后和中国国界贴合，但沿海岸线出现锯齿状空洞，小岛被切掉。原因：cutline是面矢量，栅格像元的中心点落在面外就会被裁掉；当像元一半在海一半在陆地时，crop_to_cutline按面积权重决定保留，边缘像元面积不足一半被丢弃。解决：给gdalwarp加上-wo CUTLINE_ALL_TOUCHED=TRUE，让只要与边界有接触的像元都保留；或者先裁剪到外接矩形，再做一次基于矢量的掩膜，优先保陆地像元。我的习惯是-crop_to_cutline配合-wo CUTLINE_ALL_TOUCHED=TRUE，再人工检查海岸线像元数量。

5.4 冬季NDVI高于夏季，时序曲线在1月出现尖峰

现象：黑龙江大兴安岭地区1月NDVI超过0.7，7月反而只有0.5。原因：冬季积雪在可见光波段反射率高、近红外反射率也高，NDVI公式会算出一个伪高值；如果pixel reliability没有成功标记冰雪，最大合成会把这个假高值保留。解决：在月度合成前用reliability <= 1剔出冰雪，或者在年值统计时把NDVI小于0.05的像元归为裸土/雪。注意不能全图统一设高阈值，南方常绿林冬季NDVI本来就高，阈值要分区域验证。

5.5 与站点实测FVC相关性极低，散点图一团糟

现象：用站点尺度的植被覆盖度与1km NDVI对比，R²只有0.1。原因：1km像元是一个混合像元，站点周围可能同时有树、耕地、水面和道路，站点实测的代表范围远小于1km；另外站点经纬度与像元中心存在几何偏移。解决：取站点周围7×7窗口内有效NDVI均值，而不是单点值；同时检查站点坐标是否落在水体像元上，必要时参考更高分辨率的土地覆盖数据。这个坑不是数据问题，是尺度问题，换谁跑都一样。

6. 进阶：交付前的质量验证与统计报告

6.1 用统计报告快速判断数据是否跑偏

数据做完不能直接交差，至少要出一份质量统计。我每次跑完都会生成一个JSON，把年值图的均值、标准差、最值和有效像元比例记录下来。

import json import rasterio with rasterio.open('/data/work/NDVI_2010_annual_max.tif') as src: arr = src.read(1).astype('float32') * 0.0001 valid = arr > -0.3000 report = { 'mean': float(arr[valid].mean()), 'std': float(arr[valid].std()), 'min': float(arr[valid].min()), 'max': float(arr[valid].max()), 'valid_ratio': float(valid.mean()), } with open('/data/work/NDVI_2010_raw_stats.json', 'w') as fp: json.dump(report, fp, indent=2)

逻辑说明：astype('float32') * 0.0001完成单位换算，valid掩膜排除nodata，valid_ratio表示有效像元占比。如果valid_ratio低于0.8，说明QA滤波过狠或拼接遗漏，需要回头检查。这个简单的JSON能让你在半年后回看时一眼就知道这批数据的状态，比翻文件名可靠得多。

6.2 与MOD13A1官方产品做点位交叉验证

更高一级的验证是随机抽1000个像元，与MOD13A1月度产品对应点位对比。计算双方NDVI的R²和平均绝对误差，R²大于0.9说明你的处理链没有引入明显偏差。这类验证脚本要保留，因为每次换年份或换区域都要重跑，参数只改文件路径即可。

一个长期有效的习惯：把处理参数写进输出文件名或元数据。比如NDVI_2010_annual_max_epsg4326_qa1.tif，其中qa1表示用了reliability <= 1的滤波级别。这个习惯救过我多次——同样是年度最大合成，滤波级别不同结果差异很大，文件名里不写，三个月后自己也会忘。希望帮到你。

本文还有配套的精品资源，点击获取