MODIS 2015年中国1km NDVI数据集处理与FVC计算全流程

简介本资源为2015年中国区域1km分辨率植被指数NDVI空间分布数据集面向遥感、生态、农业及地理信息领域的科研人员与学生可用于植被覆盖监测、土地退化评估、气候响应分析等场景。数据源自NASA MODIS系列MOD13A3产品经提取子数据集、拼接、投影栅格、单位换算与裁剪后生成逐月1km NDVI再以最大合成法得到年度数据投影为Albers等面积投影椭球WGS84中央经线105°标准纬线25°与47°。压缩包共5个文件约20.01MB包含1个tif栅格主文件、1个tfw坐标文件、2个xml元数据及1个txt说明文档便于直接加载与坐标校准。已有312人学习下载适合需要中国全域年度NDVI底图进行空间分析与制图的用户可快速接入ArcGIS、QGIS等平台开展后续研究。1. MODIS 2015年中国1km植被指数NDVI空间分布数据集从下载到出图的完整路径拿到一个「MODIS 2015年中国1km植被指数NDVI空间分布数据集」很多人第一反应是直接打开看数值结果发现投影不对、像元值范围诡异、边界裁切错位。这个数据集的核心价值在于它把 NASA 的 MOD13A3 月度 NDVI 产品做了中国区拼接、投影转换和最大值合成最终输出为 1km 分辨率的年际空间分布栅格。适合做全国尺度植被覆盖评估、生态质量监测、以及 NDVI 计算 FVC植被覆盖度的下游任务。但如果你不清楚它的原始投影是正弦曲线、不清楚像元值的缩放因子是 0.0001、不清楚中国边界矢量该用哪个版本后面每一步都会翻车。这篇笔记按「拿到数据 → 验证 → 裁切 → 出图 → 进阶」的顺序把每个环节的参数和坑讲透。2. 先搞懂 MOD13A3 到中国区 1km NDVI 的加工链路2.1 原始产品是什么、为什么选 MOD13A3 而不是 MOD13Q1MODIS 植被指数产品有好几个分支常见的是 MOD13Q116 天合成250m和 MOD13A3月度合成1km。这个数据集选的是 MOD13A3原因很直接做全国 2015 年一整年的空间分布如果用 250m 分辨率全国范围光单月数据量就很大拼接和后续计算对内存和磁盘都是考验。1km 分辨率下全国栅格大约 6000×4000 像元单波段浮点数据不到 100MB处理起来轻快很多。MOD13A3 本身提供的是月度 NDVI 和 EVI每个像元值经过缩放实际 NDVI DN × 0.0001。有效范围是 -2000 到 10000对应 -0.2 到 1.0。水体、云、雪等会被标记为特定填充值比如 -3000。如果你拿到数据后直接统计均值不先把填充值剔掉结果会偏得离谱。另一个关键点是投影。MOD13A3 原始文件是正弦曲线投影Sinusoidal每个瓦片覆盖 1200×1200 像元。中国区域大概涉及 h26v04、h26v05、h27v04、h27v05 等瓦片。这个数据集已经做了拼接和重投影输出通常是 WGS84 地理坐标系或 Albers 等面积投影。你需要先确认手里这份数据到底是哪种不然后面裁切中国边界时会对不上。2.2 用 Python 读取数据集并检查元信息拿到数据后别急着算先做三件事看投影、看像元值范围、看有效像元占比。下面这段代码用 rasterio 和 numpy 完成基础检查。import rasterio import numpy as np # 打开数据集注意路径换成你本地的 with rasterio.open(china_ndvi_2015_1km.tif) as src: print(CRS:, src.crs) # 确认坐标系 print(Shape:, src.shape) # 行列数 print(Bounds:, src.bounds) # 地理范围 print(Resolution:, src.res) # 像元大小 print(Nodata:, src.nodata) # 填充值 ndvi_raw src.read(1) # 读第一波段 valid ndvi_raw[ndvi_raw -2000] # 剔除填充值 print(有效像元占比:, valid.size / ndvi_raw.size) print(原始DN范围:, ndvi_raw.min(), ndvi_raw.max()) print(实际NDVI范围:, valid.min() * 0.0001, valid.max() * 0.0001)逻辑说明先读元信息确认坐标系和分辨率再用 -2000这个阈值把填充值排除。参数上缩放因子 0.0001 是 MOD13A3 的固定值不要改成 0.001 或 1。如果src.nodata返回 None说明数据集没有显式声明填充值你需要手动按 -3000 处理。如果有效像元占比低于 60%可能是中国边界裁切时把大量境外区域也算进来了或者原始瓦片拼接时没对齐。这时候要回去检查拼接步骤。2.3 中国边界矢量从哪来、怎么裁切裁切中国区 NDVI 需要一份中国国界矢量。常见做法是用国家基础地理信息中心的 1:400 万或 1:100 万标准地图或者用 GADM 的中国边界。注意 GADM 的边界包含港澳台和南海诸岛裁切时如果只想要大陆区域需要额外筛选。裁切用 rasterio 的 mask 功能最直接import geopandas as gpd from rasterio.mask import mask # 读取中国边界矢量 china gpd.read_file(china_boundary.shp) # 确保矢量坐标系和栅格一致不一致就转 with rasterio.open(china_ndvi_2015_1km.tif) as src: if china.crs ! src.crs: china china.to_crs(src.crs) # 裁切 out_image, out_transform mask(src, china.geometry, cropTrue) out_meta src.meta.copy() out_meta.update({ height: out_image.shape[1], width: out_image.shape[2], transform: out_transform }) # 写出裁切结果 with rasterio.open(china_ndvi_2015_clipped.tif, w, **out_meta) as dest: dest.write(out_image)参数说明cropTrue会把输出范围收紧到矢量边界的外接矩形减少空白像元。china.geometry是矢量几何列表如果边界是多部件要素mask 会自动处理。注意裁切后 nodata 区域会变成 0 或原 nodata 值后续统计前要再过滤一次。提示如果裁切后发现边界处像元被截断检查矢量的坐标系是否和栅格严格一致。常见错误是矢量用 CGCS2000 而栅格用 WGS84两者虽然接近但不完全等同边界处会有偏移。3. 从 NDVI 到 FVC像元二分模型的计算与验证3.1 像元二分模型的参数怎么定NDVI 计算 FVC 最常用的公式是FVC (NDVI - NDVI_soil) / (NDVI_veg - NDVI_soil)其中 NDVI_soil 是纯裸土像元的 NDVI 值NDVI_veg 是纯植被像元的 NDVI 值。这两个参数没有固定标准需要根据研究区和年份来定。常见做法是取累计频率 5% 和 95% 的分位数作为 NDVI_soil 和 NDVI_veg。import numpy as np # 假设 ndvi 是已经剔除填充值后的一维数组 ndvi_valid ndvi_raw[ndvi_raw -2000] * 0.0001 ndvi_soil np.percentile(ndvi_valid, 5) ndvi_veg np.percentile(ndvi_valid, 95) print(NDVI_soil:, ndvi_soil) print(NDVI_veg:, ndvi_veg) # 计算 FVC fvc (ndvi_valid - ndvi_soil) / (ndvi_veg - ndvi_soil) fvc np.clip(fvc, 0, 1) # 限制在 0-1 之间逻辑说明用 5% 和 95% 分位数是为了排除极端异常值的影响。如果研究区植被稀疏NDVI_veg 可能偏低这时候可以改用 90% 分位数。参数调整后 FVC 的绝对数值会变但空间分布格局通常稳定。3.2 用散点图和直方图验证 FVC 结果算完 FVC 不要直接出图先做两个验证一是看 FVC 和 NDVI 的散点关系是否单调二是看 FVC 直方图是否在 0 和 1 附近有异常堆积。import matplotlib.pyplot as plt fig, axes plt.subplots(1, 2, figsize(12, 5)) # 散点图NDVI vs FVC axes[0].scatter(ndvi_valid[::100], fvc[::100], s1, alpha0.3) axes[0].set_xlabel(NDVI) axes[0].set_ylabel(FVC) axes[0].set_title(NDVI-FVC 散点) # 直方图 axes[1].hist(fvc, bins50, edgecolorblack) axes[1].set_xlabel(FVC) axes[1].set_ylabel(像元数) axes[1].set_title(FVC 分布) plt.tight_layout() plt.savefig(fvc_validation.png, dpi150)如果散点图出现明显的分段或平台说明 NDVI_soil 和 NDVI_veg 取值不合理。如果直方图在 0 或 1 处有尖峰说明有大量像元被 clip 到边界可能需要调整分位数。3.3 和已有 FVC 产品做交叉对比如果你手头有 MODIS 的 FVC 产品比如 MOD44B 或 MOD15A2H 的 LAI/FPAR可以抽样对比。常见做法是随机选 500 个像元算两者的相关系数和 RMSE。from scipy.stats import pearsonr from sklearn.metrics import mean_squared_error # 假设 fvc_ref 是参考产品的 FVC 值 mask (fvc 0) (fvc_ref 0) r, _ pearsonr(fvc[mask], fvc_ref[mask]) rmse np.sqrt(mean_squared_error(fvc[mask], fvc_ref[mask])) print(f相关系数: {r:.3f}) print(fRMSE: {rmse:.3f})相关系数低于 0.7 时要检查两者的空间分辨率是否一致、时间窗口是否对齐。RMSE 高于 0.15 时可能是 NDVI_soil 和 NDVI_veg 的取值差异导致的系统性偏差。4. 避坑与排查NDVI 数据集处理中最容易翻车的 5 个点4.1 像元值没乘缩放因子统计结果全错现象直接读栅格得到的值在 -3000 到 10000 之间算出来的均值是几千。 原因MOD13A3 的 DN 值需要乘以 0.0001 才是真实 NDVI。 解决读取后立即做ndvi dn * 0.0001并在文档里标注缩放因子。4.2 填充值没剔除水体被算成低植被现象FVC 在湖泊、河流区域出现大面积 0 值。 原因填充值 -3000 被当作有效 NDVI 参与计算。 解决用ndvi -2000过滤或者用src.nodata指定的值过滤。注意有些数据集用 -9999 作为 nodata。4.3 投影不一致导致裁切错位现象裁切后的中国边界和 NDVI 影像对不上偏移几十公里。 原因矢量用地理坐标系栅格用投影坐标系或者两者基准面不同。 解决用gpd.read_file读矢量后先to_crs(src.crs)统一坐标系再裁切。4.4 月度合成时没做最大值合成年际 NDVI 偏低现象2015 年 NDVI 均值比相邻年份低 0.1 以上。 原因直接对 12 个月求平均把冬季低值也算进去了。 解决用最大值合成MVC取每个像元 12 个月中的最大值作为年 NDVI。如果数据集已经做了 MVC检查元数据里的合成方法说明。4.5 内存不够导致大区域处理中断现象处理全国 1km 数据时 Python 进程被 kill。 原因一次性读入整个栅格数组内存占用超过可用容量。 解决用 rasterio 的窗口读取windowed reading分块处理。或者先裁切到研究区再读入。5. 进阶技巧用 GDAL 命令行做批量重投影和格式转换如果你要处理多个年份的 MODIS NDVI 数据用 Python 脚本逐个跑太慢。GDAL 命令行工具可以批量重投影和格式转换速度比 Python 循环快很多。# 批量将正弦投影转为 WGS84并输出为 GeoTIFF for f in *.hdf; do gdalwarp -t_srs EPSG:4326 \ -tr 0.01 0.01 \ -r bilinear \ -of GTiff \ $f ${f%.hdf}_wgs84.tif done参数说明-t_srs EPSG:4326指定目标坐标系为 WGS84。-tr 0.01 0.01设置输出像元大小为 0.01 度约等于 1km。-r bilinear用双线性插值适合连续型 NDVI 数据。如果数据量特别大可以加-multi -wo NUM_THREADSALL_CPUS开启多线程。另一个技巧是用gdal_calc.py做 NDVI 缩放和填充值替换gdal_calc.py -A input_ndvi.tif \ --outfilendvi_scaled.tif \ --calcA*0.0001*(A-2000) \ --NoDataValue0 \ --typeFloat32这条命令把 DN 值乘以 0.0001同时把小于等于 -2000 的像元置为 0。--NoDataValue0指定输出文件的 nodata 值为 0方便后续在 QGIS 里设置透明。我自己的习惯是拿到任何 MODIS 数据先用gdalinfo看一眼元信息确认投影、分辨率、nodata 值再决定用 Python 还是 GDAL 命令行处理。这个习惯帮我省了很多后悔药。希望帮到你。本文还有配套的精品资源点击获取