企业官网建设流程全解析

简介：2023年中国区域一公里分辨率植被指数（NDVI）空间分布数据集，由NASA的MODIS MOD13A3月度产品经过子数据集提取、拼接、投影转换、单位换算与裁剪处理，再以最大合成法合成年度NDVI，可供生态遥感、植被覆盖监测、气候与环境变化研究者及GIS应用人员直接使用。压缩包共10个文件，包括2个TIFF栅格数据，辅以TFW空间定位文件、OVR金字塔文件、XML元数据文件和一个TXT说明文档，整体大小约67.91MB，可在常见GIS平台中快速加载显示并查看完整空间参考信息。目前已有88人学习下载，可作为区域植被动态研究的基础数据。该数据为2023年中国范围、一公里分辨率、年尺度的NDVI成果，投影采用Albers等积圆锥投影（中央经线105度，双标准纬线25度与47度，WGS84椭球），附带原始引用来源，用户无需重复下载与预处理，可直接进行空间制图、区域对比、变化检测或作为模型输入，显著节省数据准备时间。

1. 一个zip包背后的MODIS NDVI生产链路

打开压缩包，发现里面不是十几个分景HDF，而是一个已经拼好、裁好、换算过单位的china_NDVI1km_2023.tif，外加.tfw、.ovr、.xml、.aux.xml几个辅助文件。很多人拿到手直接拖进GIS，看到坐标系是Albers而不是经纬度，就开始怀疑数据坏了。实际上，这份数据是把NASA的MOD13A3月度1km NDVI产品，经过子数据集提取、跨轨道拼接、重投影、单位换算、边界裁剪，再用最大合成法（MVC）压缩成2023年年度NDVI的结果。它的价值在于省掉了每个月下载、拼接、投影的重复劳动，但也意味着你必须理解它和原始MOD13A3之间的差异：年度合成反映的是“年内最佳植被状态”，而不是全年平均值。适合需要做全国尺度生态质量评价、农业长势分析、土地利用变化研究的人，尤其是那些不想从HDF开始折腾，又担心现成数据坐标系搞错的研究生和工程师。

2. MOD13A3预处理：从月度像元到年度NDVI的取舍

2.1 先看懂MOD13A3的存储结构

MOD13A3是NASA LP DAAC发布的月度L3植被指数产品，空间分辨率1km，全球按正弦投影（Sinusoidal）分瓦片存储。每个HDF文件内部不是只有一张栅格，而是包含多个科学数据集（SDS）。直接双击HDF往往只能看到图层列表，看不到NDVI真正存在哪个子数据集里，所以第一步必须搞清楚结构。

我用gdalinfo查看HDF文件时，通常会先列出全部子数据集名称，再决定提哪一层。常见的MOD13A3 V006内部结构里，NDVI数据在“1 km monthly NDVI”这个SDS中，同时还有EVI、Pixel Reliability等。下表是几个关键子数据集：

Pixel Reliability 常被忽略，但它在做年度合成时很有用。比如某个像元在一月被云覆盖，NDVI值会异常低，而Pixel Reliability对应位会标记为“cloud”。如果你直接拿原始NDVI做最大值合成，质量差的像元也会参与竞争，得到的结果可能偏绿。稳妥的流程是把质量标记参与筛选，或者至少看一眼质量像元的空间分布。

2.2 提取、拼接与投影栅格的常见做法

原始MOD13A3是正弦投影，直接裁剪成中国范围后，如果目标坐标系不统一，后续面积统计会出大问题。数据集最终使用Albers等积投影，是因为中国国土横跨多个经度带，双标准纬线设计能让面积变形最小。常见做法是先抽子数据集，再拼接，最后用gdalwarp一次性完成投影和分辨率设置。

# 抽取某个月NDVI子数据集 gdal_translate -of GTiff \ "HDF4_EOS:EOS_GRID:\"MOD13A3.A2023001.hdf\":MODIS_Grid_1km_VI:\"1 km monthly NDVI\"" \ ndvi_202301_sin.tif # 多个瓦片拼接 gdal_merge.py -o ndvi_202301_china_sin.tif \ -n -3000 -init -3000 \ ndvi_202301_sin_tile1.tif ndvi_202301_sin_tile2.tif # 重投影到Albers，输出1000m像元 gdalwarp -t_srs "+proj=aea +lat_1=25 +lat_2=47 +lon_0=105 +datum=WGS84 +units=m" \ -tr 1000 1000 -r near -dstnodata -3000 \ ndvi_202301_china_sin.tif ndvi_202301_aea.tif

第一条命令里的HDF4_EOS:EOS_GRID是GDAL访问EOS HDF的虚拟路径，不同版本的产品内部网格名有差异，所以执行前先用gdalinfo MOD13A3.A2023001.hdf确认SDS名称。第二条命令的-n -3000表示把NoData值设为-3000，避免拼接边缘产生黑色条纹；第三条命令的-r near是最邻近重采样，对NDVI这类分类性质不太强的连续变量，一般也够用。如果希望平滑一些，可以用-r bilinear，但要注意NoData边缘会被插值污染。

2.3 单位换算和裁剪

MOD13A3的原始数值是int16，NDVI存储范围为-2000到10000左右，必须乘以0.0001才是真实NDVI。不少人拿到HDF后直接范围拉伸看图，看到数值几百以为异常，其实只是没有换算比例因子。换算和裁剪可以一步做完，也可以用gdal_calc.py单独执行：

# 裁剪中国边界并投影到最终坐标系 gdalwarp -cutline china_boundary.shp -crop_to_cutline -dstnodata -3000 \ ndvi_202301_aea.tif ndvi_202301_cn.tif # 将整数NDVI换算为浮点值 gdal_calc.py -A ndvi_202301_cn.tif \ --outfile=ndvi_202301_cn_scaled.tif \ --calc="A*0.0001" \ --NoDataValue=-0.3

这里有个细节：-3000 * 0.0001 = -0.3，所以换算后的NoData要显式写成-0.3，否则后面做最大合成时，背景值-0.3会参与计算，导致边缘出现很多“假NDVI=-0.3”的像元。裁剪边界建议使用标准Albers投影下的中国省界或国界矢量，而不是先按经纬度裁剪再投影，否则边界像元会错位。

3. 最大合成法（MVC）实现：为什么能消除云和噪声

3.1 MVC不是平均，而是选出最佳状态

月度NDVI产品虽然已经做过16天合成，但云、气溶胶、太阳高度角变化仍然会让某个月的NDVI明显偏低。最大合成法（Maximum Value Composite，MVC）对每个像元取12个月的最大值，原理是：云和气溶胶只会让NDVI变小，很难让NDVIVI变大，所以取最大值能有效保留“真实植被最旺盛”的信号。这个思路简单，但如果只看单月合成，很容易把雪、水体的高反射误判为高植被。

MVC的边界条件也很清楚：它假设年内同一个像元的真实NDVI不会出现急剧下降又回升的情况。对落叶阔叶林，夏季峰值明显，MVC结果能代表生长旺季；对常绿林，全年NDVI平稳，MVC结果也不错；但对一年两熟农田，冬小麦和夏玉米都有各自峰值，MVC只保留全年最高一次，会丢掉另一茬作物的信号。所以做农业生产时，建议同时保留各月数据，年度值只能作为补充。

3.2 用Python实现像元级最大合成

我已经把12个月的tif文件放在同一个目录，文件名格式例如ndvi_2023_01.tif到ndvi_2023_12.tif，写入前统一处理NoData：

import glob import numpy as np import rasterio monthly_files = sorted(glob.glob("ndvi_2023_*.tif")) with rasterio.open(monthly_files[0]) as src: meta = src.meta.copy() rows, cols = src.height, src.width # 用 -9999 作为内部 NoData，避免与真实NDVI混淆 stack = np.full((len(monthly_files), rows, cols), -9999, dtype=np.float32) for i, f in enumerate(monthly_files): with rasterio.open(f) as src: arr = src.read(1).astype(np.float32) arr[arr <= -0.3] = -9999 # 把背景值统一替换 stack[i] = arr # 使用掩码数组计算最大值，NoData不参与 valid_mask = stack > -9000 stack_masked = np.ma.masked_where(~valid_mask, stack) annual_max = stack_masked.max(axis=0) # 写回文件 annual_max = annual_max.filled(-9999).astype(np.float32) meta.update({"dtype": "float32", "nodata": -9999}) with rasterio.open("china_NDVI1km_2023.tif", "w", **meta) as dst: dst.write(annual_max, 1)

这段代码的核心是把NoData统一成一个极大负数，再用np.ma.masked_where把无效像元隐藏，max(axis=0)只对有效像元求最大值。如果不做掩码，直接用stack.max(axis=0)，背景值0或-0.3会把大片非植被区算成NDVI=0，最终影像看起来像蒙了一层灰。-9999不会和真实NDVI冲突，因为NDVI浮点范围几乎都在-0.2到1.0之间，最后写文件时nodata设为-9999，ArcGIS和QGIS都能正确识别。

3.3 用ArcGIS栅格计算器或arcpy实现

不写Python的可以用ArcGIS的Max工具。前提是12个月的NoData已经处理成统一的-9999或物理背景值，否则栅格计算器遇到NoData会直接返回NoData，整个中国区域都会变成空值。

import arcpy from arcpy.sa import * arcpy.env.workspace = r"D:\ndvi_2023" arcpy.env.extent = "MAXOF" arcpy.env.cellSize = 1000 pre_files = ["ndvi_202301_pre.tif", "ndvi_202302_pre.tif", "ndvi_202303_pre.tif", "ndvi_202304_pre.tif", "ndvi_202305_pre.tif", "ndvi_202306_pre.tif", "ndvi_202307_pre.tif", "ndvi_202308_pre.tif", "ndvi_202309_pre.tif", "ndvi_202310_pre.tif", "ndvi_202311_pre.tif", "ndvi_202312_pre.tif"] out_mvc = Max(pre_files) out_mvc.save("china_NDVI1km_2023.tif")

Max工具底层就是像元级最大值，但需要注意两点：第一，所有输入的像元大小和范围必须完全一致，所以前面用gdalwarp时统一-tr 1000 1000很重要；第二，如果某个像元在所有12个月都是NoData，输出也会是NoData，不会硬塞一个0进去。ArcGIS的env.extent = "MAXOF"能保证输出范围覆盖所有输入，而不会只按第一个文件的范围。

4. 文件、坐标与数据集自检：tfw、ovr、xml里有什么

4.1 解压后每个文件干什么

拿到china_NDVI1km_2023.zip，解压后会看到多个同名不同后缀的文件。很多初学者只认.tif，其他文件一律无视，这在大多数场景下没问题，但如果遇到ArcGIS读取异常、金字塔丢失或坐标系跳错，就需要知道它们的作用：

.tfw是简化版的地理定位文件，文本格式，记录像元大小、旋转量和左上角坐标。如果GeoTIFF内部坐标丢失，GIS会自动读取.tfw。.ovr是外包金字塔，能让你在缩放到全国范围时不用读取全部像元；.aux.xml里通常保存了NoData值和灰度统计，GRASS、QGIS等GDAL系程序读取时会优先使用它。删除这些文件不会损坏数据，但会拖慢首次渲染速度，也可能导致颜色渲染异常。

4.2 用代码验证投影坐标系

这个数据集的投影坐标系是Albers等积圆锥投影，椭球WGS84，中央经线105度，标准纬线25度和47度，变形比例1.0。为什么不是常见的经纬度WGS84？因为等积投影能保证面积测量准确，全国尺度统计土地覆盖面积时，用Albers比用经纬度栅格可靠得多。用projection单位是米，像元大小1000x1000，直接统计像元数再乘以1000×1000就得到面积，不用做任何换算。

import rasterio with rasterio.open("china_NDVI1km_2023.tif") as src: print("CRS:", src.crs) print("Transform:", src.transform) print("Size:", src.width, src.height) print("Bounds:", src.bounds)

我一般会先看src.crs输出是否为+proj=aea +lat_1=25 +lat_2=47 +lon_0=105 +datum=WGS84，再手动确认transform里的像元尺寸接近1000,1000。如果看到+proj=longlat或degree单位，说明投影信息丢失或被人为覆盖。这时可以手动在GIS里把坐标系定义为上述Albers参数，但要注意lat_0通常默认0，标准纬线顺序不影响结果。

4.3 容易踩的坑：NoData和值域

最大合成后的NDVI理论上应该在-0.2到1.0之间，但水体、冰雪、裸地有时会出现很低的负值，甚至-0.3。如果你在ArcGIS里看到整个图层都是黑色，先打开符号系统，看拉伸范围是否把NoData和真实值混杂在一起。.aux.xml里如果有<MDI key="NODATA">-9999</MDI>，一部分GIS版本会正确识别，另一部分版本需要你在图层属性里手动设置NoData为-9999。

另外要注意：.tif内部有内嵌坐标，.tfw只是备用，如果你把.tif拷贝到新的文件夹却忘了同时拷贝.tfw，绝大多数情况不影响使用，但当你用某些老旧的C++程序直接按TFW读文件时，就可能出现仿射参数错位。所以日常交换数据时最好保持zip包完整，不要只挑选.tif发送。

5. 快速验证年度NDVI数据的三个实用技巧

5.1 用Python看分布并做合理性质检

拿到年度NDVI后，第一步不是画图，而是计算有效值范围和均值。如果均值在0.3到0.7之间，说明中国大部分区域植被覆盖较好；如果均值接近0.8，大概率有异常高值参与了合成。

import rasterio with rasterio.open("china_NDVI1km_2023.tif") as src: ndvi = src.read(1, masked=True) print("Valid sample count:", ndvi.count()) print("Min:", float(ndvi.min()), "Max:", float(ndvi.max())) print("Mean:", float(ndvi.mean()))

masked=True会自动把NoData排除在统计外，避免背景值拉低均值。合理的年度NDVI最大值不应超过1.0，但因为有水体镜面反射或积雪，可能略超，如果看到2.0这样的数值，基本可以断定NoData没设对。

5.2 用pyproj把经纬度转成Albers坐标后采样

很多场景需要比较“某个县”或“某个野外台站”的NDVI值。台站经纬度是WGS84，但栅格坐标是Albers，直接读取很难得到准确结果。用pyproj转换后采样最稳：

from pyproj import Transformer import rasterio transformer = Transformer.from_crs("EPSG:4326", "+proj=aea +lat_1=25 +lat_2=47 +lon_0=105 +datum=WGS84", always_xy=True) lng, lat = 108.36, 34.21 x, y = transformer.transform(lng, lat) with rasterio.open("china_NDVI1km_2023.tif") as src: for val in src.sample([(x, y)]): print("NDVI at Qinshan station:", val)

always_xy=True保证输入输出都是经度在前、纬度在后，不然transform函数默认按纬度经度顺序，容易得到一组颠倒且不报错的坐标。采样输出是单值数组，如果返回的是-9999，说明你采样的点正好落在NoData像元，需要扩大邻域或检查坐标范围是否在国界内。

5.3 与往年数据比较时先统一投影和分辨率

年度NDVI最适合做年际变化检测，但前提是两年数据必须在同一投影、同一像元大小、同一覆盖范围下比较。不要拿2023年的Albers数据和2020年的经纬度坐标数据直接相减，结果是位移整公里量级的噪声。

gdalinfo -stats china_NDVI1km_2023.tif

通过gdalinfo -stats检查输出中的STATISTICS_MINIMUM、STATISTICS_MAXIMUM和STATISTICS_MEAN，如果均值在两个相邻年份间出现超过0.2的跳跃，不要急着归因于气候，先检查源数据是不是存在传感器退化或质量标记问题。验证完成后，再把china_NDVI1km_2023.tif重命名归档，保留原始投影参数和NoData设置，后续任何分析都从这个基准出发。

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