企业官网建设流程全解析

简介：这是一份基于NASA MOD13A3产品加工而成的2024年中国1km植被指数（NDVI）空间分布数据集，适合生态学、地理学、农学等领域的师生和科研人员用于植被覆盖监测、土地利用变化分析及遥感制图。数据采用Albers等积圆锥投影，空间分辨率1km，时间覆盖2024年全年，经过子数据集提取、拼接、投影、裁剪和最大合成等步骤生成，可直接在ArcGIS、QGIS等软件中加载使用。资源包共10个文件，以GeoTIFF栅格数据为主，辅以tfw坐标信息、xml元数据和ovr金字塔文件，压缩包大小约75.87MB，便于快速下载与处理。目前已有361人学习下载。借助该数据集，用户可以快速获取全国尺度年度NDVI空间分布，免去原始HDF数据批量下载与预处理流程，直接用于植被长势评估、干旱监测、生态环境评价等研究或教学展示。

1. MODIS 2024年中国1km NDVI：一台机器、一天下载、一周处理的地表植被底图

做生态评估、农业长势监测或者城市绿化分析的人，2024年快结束的时候大概率会撞上同一个需求：要一份全国范围、1km分辨率、逐月更新的植被指数（NDVI）数据。MODIS系列的NDVI产品是目前唯一能免费满足这个条件的长时序数据源，但真正从NASA拉过数据的人都知道，所谓“2024年中国1km NDVI数据集”不是一个能直接下载的现成文件，而是要自己从十几个分幅（tile）的HDF文件里拼出来的一套数据。这篇笔记就是把这条路径完整讲清楚：选什么产品、怎么下载、怎么把HDF转成能用的GeoTIFF、质量控制怎么做、坑在哪，最后怎么把NDVI进一步算成植被覆盖度（FVC）。读者如果是做过遥感但没碰过MODIS的人，跟着做一遍就能拿到可用的全国逐月NDVI；如果已经拉过数据但被投影或者质量波段卡住，直接跳到第5章排查。

2. 选MOD13A3还是MOD13Q1：月度1km产品的参数权衡与下载前准备

2.1 两个候选产品在时间分辨率上的根本差异

MODIS植被指数产品家族里有几个经常被搞混的名字：MOD13Q1是250m、16天合成的NDVI/EVI，MOD13A1是500m、16天合成，MOD13A2是1km、16天合成，MOD13A3是1km、月度合成，MYD13A3是Aqua卫星对应的月度版本。标题里写的是“1km植被指数空间分布数据集”，而且按年份来看，最省事的做法是选MOD13A3，因为一年12个月、每个月一个文件，处理量小，时间分辨率对齐了日历月，拿来做年度趋势、季节曲线这类分析非常顺手。MOD13Q1虽然空间分辨率高得多，但一年23期（16天合成加年末一期），全国要下载二十多个tile乘以23期，处理时间翻好几倍，而且250m的产品在实际生态评估中并不一定带来更高的精度——合成窗口短意味着云污染残留更多，质量控制工作量反而更大。

从算法上看，MOD13A3不是独立反演的，它是在MOD13A2（1km、16天）的基础上做的月度合成，具体做法是把一个月内所有可用观测做加权平均，权重来自合成期的质量标记。因此MOD13A3的NDVI值本身已经带了一定程度的去云和去气溶胶处理，但这不代表可以直接拿来用，后面第4章会细说质量波段的事。还有一点容易被忽略：MOD13A3的NDVI是Terra卫星的，MOD13A2是Terra的16天产品，MYD13A3是Aqua的月度产品。Terra上午过境、Aqua下午过境，两者对同一块地的观测角度和云干扰不同，直接对比的时候要注意来源一致，混用会在时间序列上造出额外噪声。

2.2 LP DAAC的tile覆盖：中国区域需要哪几个分幅

MODIS采用正弦投影（Sinusoidal），全球按10度乘10度划分成h（水平）和v（垂直）编号的分幅。中国地跨约73°E到135°E、18°N到53°N，对应到分幅上，常见覆盖范围大致是h23v04、h24v03到h28v06这一片。我处理2024年数据时用到的分幅集合是：h23v04，h23v05，h24v03，h24v04，h24v05，h25v03，h25v04，h25v05，h25v06，h26v03，h26v04，h26v05，h26v06，h27v04，h27v05，h27v06，h28v05，h28v06，共18个分幅。实际项目里建议不要凭记忆定列表，而是先把中国国界转成正弦投影坐标系，算出边界框，再去NASA的Earthdata检索页面勾选tile，这样能避免漏掉海岸线附近的小分幅。

下载之前要做的准备工作有三件。第一，注册NASA Earthdata账号并创建应用令牌，因为批量下载脚本依赖这个令牌做HTTP鉴权，手动在浏览器里一个个点下载在120个月度文件面前完全不现实。第二，明确数据版本，目前MODIS陆地产品主流是Collection 6.1，文件命名里会带061字段，选这个版本，不要用老的C6。第三，想清楚要Terra还是Aqua，或者两个都要。做年度分析我一般只要MOD13A3，但如果研究区云多，可以把MYD13A3也下下来做平均或者做交叉验证。

2.3 下载命令与本地目录组织的建议

用Python的earthaccess库可以省掉很多HTTP细节，下面这段代码是我常用的下载流程：

import earthaccess # 登录，会弹出交互式输入用户名和密码，或者用环境变量 earthaccess.login() results = earthaccess.search_data( short_name="MOD13A3", version="061", bounding_box=(73, 18, 135, 53), # (西, 南, 东, 北) temporal=("2024-01-01", "2024-12-31") ) # 这一步会返回数据文件链接列表，按需过滤 files = earthaccess.get_data(results, "G:/MODIS_NDVI_2024")

bounding_box参数用的是WGS84经纬度坐标，检索系统会返回覆盖这个框的所有分幅文件。temporal参数指定2024年全年，每个分幅每月一个文件，一共18×12=216个HDF文件。earthaccess.get_data会逐个下载，速度取决于网络，国内环境下每个文件大约几十MB，全下完约10GB，一晚上能跑完。

下载完的目录结构我强烈建议立即按“分幅→日期”重新组织，别让文件堆在一个目录里。文件名自带日期字段，比如MOD13A3.A2024001.h26v04.061.2024031.hdf里的A2024001表示数据观测起始日期（儒略日001），但很多人在后期排序时不解析这个字段，直接用文件名排序，结果h26v04和h26v05之间互相穿插，2024年12个月全错位。我一般会写一个小脚本，解析A2024后接的三位儒略日，转成标准化日期写到文件名前缀，再按tile分类到子目录。这个习惯能省掉后面所有时间序列处理的麻烦。

3. 用GDAL/Rasterio把HDF转成可用的GeoTIFF：拼接、投影与裁剪

3.1 HDF4的子数据集结构与读取路径

MOD13A3的HDF文件是HDF4格式，里面不是一个单独的NDVI栅格，而是多个科学数据集（Scientific Dataset，简称SDS）打包在一起。用GDAL打开时要通过子数据集路径定位，常见的子数据集名称是“1 km monthly NDVI”。第一次接触的人最容易在这里翻车：直接gdal.Open(文件名)打开的是一个HDF容器，不是栅格本身，必须按HDF4_EOS:EOS_GRID:文件名:MOD_Grid_monthly_1km_VI:1 km monthly NDVI这种格式的路径去访问。

用rasterio读取的代码如下：

import rasterio hdf_path = "MOD13A3.A2024001.h26v04.061.2024031.hdf" subdataset = f"HDF4_EOS:EOS_GRID:{hdf_path}:MOD_Grid_monthly_1km_VI:1 km monthly NDVI" with rasterio.open(subdataset) as src: ndvi = src.read(1) # 第一个波段 profile = src.profile print(profile["crs"], profile["height"], profile["width"]) # crs会显示为正弦投影(Sinusoidal)，行列数约1200x1200

这段代码的关键是subdataset字符串的拼法，不同MODIS产品内部的网格名不同，MOD13A3是MOD_Grid_monthly_1km_VI，MOD13Q1则是MOD_Grid_250m_16_days_VI，如果直接套用会报does not exist in the file system错误。打开后profile["crs"]能看到数据自带的正弦投影坐标系，注意这个投影的像元不是正南正北的矩形，直接按经纬度去裁剪一定会出问题。

顺带说一个容易忽略的事：MOD13A3里的NDVI是Int16整型存储的，真实值等于DN值乘以0.0001的比例因子。有效范围一般是-2000到10000，对应真实NDVI的-0.2到1.0，填充值（fill value）是-3000，但实际文件里也可能出现-28672之类的边界值。读取后第一件事就是把DN转成浮点真实值，并且把填充值替换成NaN，否则后面统计全是错的。

3.2 分幅拼接与正弦投影转WGS84

18个分幅覆盖中国，但分幅之间既有重叠又有缝隙，而且正弦投影的图幅在经纬度上不是规则的。处理顺序有讲究：先重投影再拼接，而不是先拼接再投影。因为如果先把18个分幅拼成一个大的正弦投影影像再转投影，中间步骤的内存占用会很高，而且边缘分幅的无效值区域会在拼接时混进来。我一般分两步走：先用gdalwarp把每个分幅单独转成WGS84，再用gdal_merge.py拼接。

命令行方式最简单，适合脚本化：

gdalwarp -t_srs EPSG:4326 -r near -dstnodata -9999 \ "HDF4_EOS:EOS_GRID:MOD13A3.A2024001.h26v04.061.2024031.hdf:MOD_Grid_monthly_1km_VI:1 km monthly NDVI" \ h26v04_2024001_wgs84.tif

这里-t_srs EPSG:4326把正弦投影转成经纬度，-r near用最近邻重采样保持原始像元值不插值，-dstnodata -9999把无效值统一设成-9999。注意不要用bilinear（双线性）重采样做NDVI，因为双线性会混合有效值和填充值，在海岸线、云边界处造出介于两者之间的伪值，后续阈值分割全部受污染。重采样方法这个细节，是很多教程不会提但实际影响很大的参数。

重投影之后的影像分辨率不再是严格的1km，而是取决于目标坐标系下的像元尺寸。转成WGS84后纬度方向的像元在赤道是0.008983度约等于1km，但在中国纬度（北纬18到53度）经度方向的物理距离会被压缩，所以单看度数分辨率没有意义。如果后续要做面积统计，我更推荐把目标坐标系设成Albers等积投影（EPSG:102025或按区域选取中央经线），这样像元面积才一致。如果只是做趋势分析或者和气象站点数据匹配，WGS84够用了。

3.3 按国家边界裁剪并统一nodata

拼接完的全国影像范围是一个覆盖中国及周边国家的大矩形，裁剪到国界边界能显著减少无关像元，后续处理速度更快，统计也更干净。裁剪用gdalwarp加-cutline参数最直接：

gdalwarp -cutline china_2024.shp -crop_to_cutline \ -dstnodata -9999 -of GTiff merged_2024001_wgs84.tif \ china_ndvi_2024001.tif

注意china_2024.shp必须是WGS84坐标，和影像保持一致。裁剪后要复查一下nodata分布，gdalwarp在裁剪边界外填入-dstnodata设定的值，但影像内部原有的填充值还是原来的-3000或者-28672，没有统一替换。这时候需要一次栅格运算把所有非有效值统一成-9999：

import rasterio import numpy as np with rasterio.open("china_ndvi_2024001.tif") as src: data = src.read(1).astype(np.float32) profile = src.profile data[data < -5000] = -9999 # 把填充值、边界外值统一 data[data < -2000] = -9999 # 低于有效范围下限的也排除 ndvi_real = data * 0.0001 # DN转真实NDVI值 ndvi_real[ndvi_real < -0.5] = np.nan # NDVI低于-0.5基本不可信，设NaN profile.update(dtype=np.float32, nodata=None) with rasterio.open("china_ndvi_2024001_clean.tif", "w", **profile) as dst: dst.write(ndvi_real, 1)

这一步的核心是执行了两次过滤：第一次清理物理存储上的填充值，第二次把物理值域映射到NDVI真实范围。这里有个经验值：NDVI低于-0.5的像元即使质量标记没问题，也大概率是水体或者沙漠中的干扰像元，在做陆地植被分析时直接设NaN比留着更安全。水体在红波段吸收强、近红外反射低，NDVI通常为负，但如果研究目标包括湿地水体，需要单独保留，不要一刀切。

4. 像元可靠性与时间序列合成：从NDVI到FVC的关键一步

4.1 pixel reliability波段：云、雪、气溶胶如何污染NDVI

MOD13A3文件里除了NDVI本身，还带了一个“pixel reliability”子数据集，这是质量控制的入口。这个波段每个像元取值从0到4，含义是：0表示高质量可靠数据（good），1表示边缘数据（marginal），2表示雪或冰覆盖，3表示云覆盖，4表示被气溶胶或阴影干扰。很多人的错误是只读NDVI波段不读质量波段，拿到的月度值里面混着大量云和雪的观测，做出来的时间曲线跳来跳去，还以为是植被的真实波动。

处理逻辑很简单：在读取NDVI的同时读取pixel reliability，把质量值大于等于2的像元全部屏蔽掉。不过要注意，月度合成产品在合成时已经做过一次优选，剩下的质量标记为1的像元不代表不可用，在云多的月份（比如华南的春季）可选的可靠像元太少，保留marginal像元能避免大片空洞。实际项目中我采用双阈值：做年度最大合成时只接受0和1，做逐月趋势分析时0和1都保留，但把1单独标记出来做敏感性分析。

屏蔽的代码实现：

with rasterio.open(subdataset_rel) as rel_src: reliability = rel_src.read(1) valid_mask = (reliability == 0) | (reliability == 1) ndvi_clean = np.where(valid_mask, ndvi_real, np.nan)

4.2 MAXVALUE合成与Savitzky-Golay滤波

月度NDVI数据虽然已经是合成产品，但受残留云和大气噪声影响，逐月曲线上经常可以看到单月骤降随后回升的“毛刺”。如果是做年度最大NDVI（简单理解为一年中植被最茂盛时期的NDVI），直接取12个月的max就行，这一步用最大值合成天然能滤掉大部分云噪声，因为云只会让NDVI变低，不会变高。

但做季节曲线或物候分析时，最大值合成帮不了忙，需要滤波平滑。我常用的方法是Savitzky-Golay滤波，它能在保留植被生长曲线峰值形态的同时去掉短周期噪声，比滑动平均好得多，滑动平均会把春季NDVI快速上升的拐点抹平。

from scipy.signal import savgol_filter # month_series: 12个月NDVI的数组，形状为 (12, height, width) # 对每个像元做时间维平滑 window_length = 7 # 必须为奇数，月度序列中7约覆盖一个生长季 polyorder = 2 # 多项式阶数，2阶足够，高阶容易过拟合噪声 smoothed = np.full_like(month_series, np.nan, dtype=np.float32) for i in range(month_series.shape[1]): for j in range(month_series.shape[2]): ts = month_series[:, i, j] valid_idx = np.isfinite(ts) if valid_idx.sum() >= 5: # 至少5个月有效，否则不滤波 filled = np.interp( np.arange(12), np.nonzero(valid_idx)[0], ts[valid_idx] ) # 先用线性插值填补缺测 smoothed[:, i, j] = savgol_filter(filled, window_length, polyorder)

这段注释里的要点是：savgol_filter不接受NaN，必须预填。缺测月份用线性插值填补，但只用于滤波前的填充，滤波后的平滑值才是输出——这相当于给插值结果加了一个先验约束，不让它自由波动。window_length=7对应7个月窗口，对于生长季大约5到8个月的植被来说，既能平滑掉单月毛刺又不会削掉峰值。北方落叶阔叶林春季NDVI从0.3涨到0.8可能只要两个月，窗口再大会把这一过程拉平，物候始期估算就会偏晚。

4.3 用NDVI分位数估算FVC植被覆盖度

NDVI的绝对数值受土壤背景、大气残留和传感器定标影响，跨区域直接比较的意义有限，所以业务上更常把NDVI转换成植被覆盖度（FVC）。这也是很多做生态评估的项目最终交付的指标。FVC的计算方法是像元二分模型：

FVC = (NDVI - NDVI_soil) / (NDVI_veg - NDVI_soil)

其中NDVI_soil是纯裸土像元的NDVI，NDVI_veg是纯植被像元的NDVI。理论上这两个端元应该从实测光谱获取，但实际项目里最可靠的获取方式是从影像自身统计——取NDVI累积分布的5%和95%分位。用百分位数而不是全局最大最小值是这里最常见的坑，全局最大最小值会被异常高值（比如云边缘、耀斑区）和异常低值（深水、阴影）拉偏，导致FVC整体压缩。用百分位数做两端截断后，把FVC限幅到0到1之间：

ndvi_5 = np.nanpercentile(ndvi_annual_max, 5) # 裸土端元 ndvi_95 = np.nanpercentile(ndvi_annual_max, 95) # 全植被端元 fvc = (ndvi_annual_max - ndvi_5) / (ndvi_95 - ndvi_5) fvc = np.clip(fvc, 0, 1)

这里用的ndvi_annual_max是第4.2节得到的12个月最大值合成结果，而不是逐月NDVI。原因是FVC估算针对的是年度最茂盛时期的覆盖度，逐月FVC受物候影响太大，年初算出来覆盖度低不是植被退化而是还没长出来。如果项目需要逐月FVC，那就要对每个月的NDVI分别做分位数提取，而不是用同一个年度端元。

5. 做这套数据最容易翻车的5个坑：现象、原因与排查顺序

5.1 坑一：HDF自带的正弦投影让经纬度裁剪完全失效

现象：按经纬度边界裁剪后得到的是全黑或只有零星像元的影像，拉伸显示完全看不出中国轮廓；或者用GeoJSON做mask时报坐标范围不重叠的错误。

原因：MOD13A3的HDF内部是正弦投影，像元在地球表面是等面积分布的，但投影后的行列方向和经纬度不是简单对应。直接拿经纬度坐标去索引像元行列号，边界框算出来的行列区间全部落在无数据的角落。

解决：先重投影到EPSG:4326再裁剪，这条在第3.2节已经强调过。补充一个排查手段：打开gdalinfo看Corner Coordinates，如果显示的单位是米而不是度，就先做gdalwarp -t_srs EPSG:4326。以后凡是拿到任何MODIS的HDF数据，第一步就是查投影，不要假设任何产品的坐标系。

5.2 坑二：把MOD13C2或MYD13A3当成1km月度数据

现象：处理完发现影像的行列数不是预期的一千多乘一千多，而是700多行或者3000多行，或者两个月份的影像边界无法对齐。

原因：MODIS植被指数产品家族里有几个容易混淆：MOD13C2是CMG（Climate Modeling Grid）产品，空间分辨率是0.05度，约5.5km，行列数只有3600×7200；MYD13A3是Aqua卫星的月度1km数据，和MOD13A3同分辨率但过境时间不同。下载检索时如果没注意平台和产品代码，很容易混入不同来源。

解决：下载后统一用gdalinfo检查行列数和crs，按MOD13A3、MYD13A3、MOD13C2分类归档。如果混用了Aqua和Terra数据做合并分析，时间序列上会出现系统性偏差（Terra上午过境，Aqua下午过境），这种偏差比随机噪声更危险，因为它有方向性，比如下下午云多时Aqua观测频率下降，NDVI系统性偏低。

5.3 坑三：冬季NDVI出现负值或伪高值，被当成植被真实信号

现象：12月和1月的NDVI在全国影像上出现大片负值（尤其华南和西南多云区），或者青藏高原部分像元出现高于夏季的反常高值，月度时间曲线剧烈抖动。

原因：这是云、雪、BRDF三重因素叠加的结果。云层在红波段反射高，会让NDVI往负值方向拉；积雪在可见光强反射，近红外吸收，计算出来也是低值；而高原地区太阳高度角低、传感器视角大，双向反射分布函数（BRDF）效应会放大某些像元的反射率信号，导致NDVI虚高。月度合成对云的抵抗力有限，尤其是冬季云多的月份。

解决：先用pixel reliability波段把雪和云（值2和3）屏蔽掉，替换成NaN。如果缺测太多，用第4.2节的线性插值补缺但不输出插值结果，只输出滤波平滑值。还有一个更彻底的办法：对于关键的分析场景，改用MCD43A4（BRDF校正后的反射率产品）自己算NDVI，它的质量更高但处理量陡增，适合研究区域不大的项目。

5.4 坑四：文件按tile编号排序导致时间序列错位

现象：逐月NDVI提取的时间序列在某几个月份突然跳到负数或者突变到异常高的值，检查单月影像发现该月份被替换成了另一个tile的编号。

原因：文件命名里既有tile编号（如h26v04）又有儒略日（如A2024289），用普通字符串排序时，文件先按tile编号排列再按日期排列，如果脚本里用glob.glob加sorted()直接处理后，12个月中每一个月的序列实际上来自同一个tile的12个不同观测，序列错位了。

解决：统一从文件名中解析儒略日字段并转成整数，按它排序。正确做法：

import re, glob def julian_sort_key(path): fname = path.split("/")[-1] match = re.search(r"A(\d{4})(\d{3})", fname) return int(match.group(2)) # 儒略日作为排序键 files = sorted(glob.glob("G:/MODIS_NDVI_2024/*/*.tif"), key=julian_sort_key)

这个坑排查起来最耗时的原因是现象和tile模糊不清，单月看数据很正常，一拉时间序列就崩。处理完排序后，建议在代码里加入检查：连续两个文件名的儒略日差值应该在28到31天之间，超出这个范围直接抛异常终止，而不是带病跑完全部12个月。

5.5 坑五：FVC估算用全局最大最小值导致覆盖度截断

现象：FVC影像上城市建成区和裸土像元的覆盖度不是接近0而是0.2甚至0.3，茂密森林像元的覆盖度到不了1，影像整体对比度很差，直方图集中在0.4到0.6之间。

原因：直接用np.nanmax和np.nanmin取NDVI端元时，全局最大值被个别异常像元（云边缘、传感器饱和）拉高，全局最小值被深水体和阴影拉低，导致FVC的线性拉伸范围比真实动态范围宽，中间值被压缩。

解决：改用分位数端元，按第4.3节的做法，用5%和95%分位数代替最大最小值。另外分位数计算时要按年度最大NDVI来做，不要对每个月分别算，否则每个月的土壤和植被端元都不一样，FVC在时间上不可比。这个坑一旦踩进去，后面做的所有生态评估指标都会带着系统性偏差，而且很难从最终结果上看出来，只能在FVC直方图上发现异常。

6. 验证与进阶：把2024年NDVI序列用起来

6.1 与MCD12Q2物候产品做交叉验证

拿到一套干净的逐月NDVI之后，第一件事不要急着画图或者算统计量，而是做交叉验证。一个低成本高收益的做法是把NDVI时间序列和MODIS的MCD12Q2物候产品对比——MCD12Q2直接给出了每个像元的生长季开始日期（SOS）和结束日期（EOS），从NDVI曲线上用阈值法（比如NDVI达到当年振幅的20%时记为返青期）提取的物候日期，理论上应该和MCD12Q2给出的日期高度相关。如果两者差异超过20天，大概率是滤波参数有问题或者质量屏蔽过度，这时候回头调window_length或者检查缺测填充逻辑。我自己的经验是，生长季中期（6到8月）的NDVI值和物候日期相关系数能到0.8以上才算数据干净。

6.2 空间统计与异常年份筛查的小技巧

2024年单年的NDVI要想看出价值，最好和过去几年的基期做对比。我没有让你非得拉十年数据，但至少可以把2024年的年度最大NDVI和2020到2023年的平均值做差，得到一张NDVI距平图。距平图比绝对NDVI图更能说明问题，因为绝对NDVI受植被类型和气候背景控制，而距平图直接显示“今年哪里比往年绿了、哪里黄了”。

具体操作上有一个小技巧：计算生长季累积FVC比单看最大FVC更稳定。把4到10月每月的FVC加总，得到一个生长季累积覆盖度，它对单月的异常噪声更不敏感，而且和作物产量、草地生产力的相关性通常更好。我用这个指标做森林恢复监测时比单月NDVI的变化检测稳定得多。2024年如果要做生态质量评价，我强烈建议输出一张生长季累积FVC而不是12个月逐月堆叠的NDVI。

最后说一个操作习惯：每一版处理好的数据都带上处理日期后缀，不要覆盖上一版。做MODIS数据处理的变量太多了——tile列表、云阈值、滤波窗口、端元分位数，每一个参数微调都会得到不同的结果，没有版本管理等于没有后悔药。我过去在季度项目里因为覆盖了上一版结果，不得不重新跑三天预处理，从那以后所有中间产物都带_v1、_v2后缀，参数记录写在同目录下的processing_log.txt里。数据做完、图出来、参数写清楚，这条链路才算真正闭环。希望帮到你。

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