简介本资源为2015年中国区域1km分辨率植被指数NDVI空间分布数据集面向遥感、生态、农业及地理信息领域的科研人员与学生可用于植被覆盖监测、土地退化评估、气候变化分析等场景。数据源自NASA MODIS系列MOD13A3产品经子数据集提取、拼接、投影转换、单位换算与裁剪后采用最大合成法生成年度NDVI投影为Albers等面积圆锥投影椭球WGS84中央经线105°标准纬线25°与47°。压缩包共5个文件约20.01MB包含1个tif栅格主文件、1个tfw坐标文件、2个xml元数据及1个txt说明文档便于直接加载与坐标配准。已有312人学习下载。读者可获取完整年度NDVI栅格用于空间分布制图、时序对比与模型输入并借助说明文档快速理解数据来源与处理流程。1. 从一份 2015 年中国 1km NDVI 数据说起它到底能干什么如果你手上正好有一份 MODIS 2015 年中国 1km 植被指数NDVI空间分布数据集第一反应大概率不是「这数据真漂亮」而是「我拿它能算什么」。我最早接触这类数据是为了做县域尺度的植被覆盖度FVC反演当时手头只有行政边界和一堆气象站点最缺的就是一张能覆盖全国、时间对得上、分辨率又不至于太粗的底图。NDVI 正好补上了这个缺口——它用近红外和红光波段的反射率差值归一化直接反映植被绿度1km 像元意味着全国大约 960 万平方公里被切成不到一千万个格子单省跑一遍内存扛得住做趋势分析也不至于被云污染撕得七零八落。这份数据集的核心价值在于「空间分布」四个字。它不是单点观测不是站点插值而是栅格化的面状覆盖能直接和土地利用、气候区划、DEM 做叠加。适合谁用做生态评估的、算 FVC 的、盯退耕还林效果的、写区域规划报告的甚至只是想给论文配一张全国植被长势图的都能从这份 2015 年的切片里挖出东西。但前提是你得先搞清楚它的投影、填充值、有效范围否则后面每一步都是坑。2. MODIS NDVI 的物理含义与 2015 年 1km 产品的选型逻辑2.1 NDVI 为什么能代表植被从波段反射到归一化差值NDVI 的公式简单到一行就能写完(NIR - Red) / (NIR Red)。但简单不代表随便用。健康植被的叶肉组织强烈反射近红外同时叶绿素吸收红光所以 NIR 高、Red 低差值大NDVI 就高。裸土、水体、云、雪则相反NDVI 趋近于零甚至为负。MODIS 传感器在 2015 年已经稳定运行了十几年波段响应函数成熟辐射定标精度在 2% 以内这是它比很多国产传感器更早被拿来做长时序分析的原因。但要注意NDVI 不是万能的。它饱和得早——当叶面积指数超过 4 左右NDVI 就涨不动了所以密林区容易低估变化。另外它受土壤背景影响大稀疏植被区信号弱。2015 年这个时间点选得也巧既避开了 2000 年代初传感器刚上天时的数据波动又没到后来部分波段衰减明显的阶段做基准年比较合适。2.2 为什么是 1km 而不是 250m 或 500m分辨率与全国覆盖的权衡MODIS 的 NDVI 产品有 250m、500m、1km 三档。250m 最细但全国跑一遍数据量直接翻四倍而且 250m 只有红光和近红外两个波段云掩膜和大气校正的辅助信息少质量控制字段不如 1km 丰富。500m 是折中但很多做省级以下单元分析的人会发现 500m 在县级尺度上还是太粗一个县没几个像元统计意义弱。1km 产品的优势在于第一它带了详细的 QA 波段能逐像元判断云、雪、阴影、气溶胶影响第二全国范围单年数据量在几百 MB 到 1GB 量级普通笔记本能处理第三很多已有的土地利用、气象、土壤数据集本身就是 1km 或更粗匹配起来不用重采样。我一般会建议做全国或大区域趋势分析1km 足够做地块级或城市内部绿化评估直接上 250m 或更高分辨率别在 1km 上硬抠。2.3 2015 年中国区数据的获取与格式确认这份数据集常见格式是 GeoTIFF 或 HDF投影多为正弦投影Sinusoidal也有已经转成 Albers 等积投影的版本。拿到手第一件事不是打开看而是用gdalinfo把元数据读出来。下面这段命令是我每次拿到新栅格必跑的gdalinfo MODIS_NDVI_2015_China_1km.tif重点看几个字段Size is后面的行列数Coordinate System is里的投影和基准面NoData Value填充值Band 1的Type和STATISTICS。如果 NoData 是 -3000 或 32767 这类值后面计算前必须掩掉否则全国均值会被拉低一大截。投影如果是正弦的做面积统计前先重投影到等积投影不然高纬度地区面积会失真。提示不要直接用 QGIS 打开就目测。先跑gdalinfo把 NoData、投影、行列数记在脚本注释里后面每一步都引用这些参数。3. 用 Python 把 2015 年 NDVI 栅格跑通从读取到 FVC 反演3.1 读取栅格与有效像元掩膜rasterio 最小代码读取用rasterio最稳它比 GDAL 的 Python 绑定更友好。下面这段代码做了三件事读数据、读 NoData、生成有效掩膜。import rasterio import numpy as np path MODIS_NDVI_2015_China_1km.tif with rasterio.open(path) as src: ndvi src.read(1).astype(float32) # 读第一波段 nodata src.nodata # 读取 NoData 值 profile src.profile # 保存元数据供后续写出 print(行列数:, src.height, src.width) print(投影:, src.crs) print(NoData:, nodata) # 构建有效掩膜排除 NoData 和 NDVI 异常值 valid np.ones(ndvi.shape, dtypebool) if nodata is not None: valid (ndvi ! nodata) valid (ndvi -1.0) (ndvi 1.0) # NDVI 理论范围 ndvi_valid np.where(valid, ndvi, np.nan) print(有效像元占比: {:.2f}%.format(valid.mean() * 100))逻辑说明src.read(1)读的是第一波段如果文件里还有 QA 波段波段顺序要确认。nodata可能是 None所以先判断。NDVI 理论范围是 -1 到 1超出这个范围的像元基本是填充或异常直接掩掉。最后用np.nan替换无效值方便后面用np.nanmean统计。参数说明astype(float32)是为了省内存全国 1km 栅格大概 8000×6000 左右float32 占不到 200MB。如果内存紧张可以分块读但 2015 年单年数据一般不用。3.2 从 NDVI 到植被覆盖度 FVC像元二分模型参数怎么定FVC 的像元二分模型公式是FVC (NDVI - NDVI_soil) / (NDVI_veg - NDVI_soil)。关键在于 NDVI_soil 和 NDVI_veg 怎么取。常见做法是取累计频率的 5% 和 95% 分位数而不是固定值。下面代码演示# 基于有效像元计算分位数 ndvi_flat ndvi_valid[~np.isnan(ndvi_valid)] ndvi_soil np.percentile(ndvi_flat, 5) # 裸土端元 ndvi_veg np.percentile(ndvi_flat, 95) # 纯植被端元 print(NDVI_soil {:.4f}, NDVI_veg {:.4f}.format(ndvi_soil, ndvi_veg)) # 像元二分模型计算 FVC fvc (ndvi_valid - ndvi_soil) / (ndvi_veg - ndvi_soil) fvc np.clip(fvc, 0, 1) # 截断到 0-1 # 统计全国均值 print(2015 年全国 FVC 均值: {:.4f}.format(np.nanmean(fvc)))逻辑说明分位数取 5% 和 95% 是经验值目的是排除极端异常像元。如果研究区是纯森林NDVI_veg 会偏高如果是干旱区NDVI_soil 可能接近 0。所以最好分区域定参数不要全国一刀切。np.clip把超出 0-1 的值截断避免负 FVC 或大于 1 的荒谬值。参数说明分位数可以调比如 2% 和 98% 更激进5% 和 95% 更保守。我一般先跑一遍看直方图如果 NDVI 分布双峰明显分位数法很稳如果单峰说明研究区植被类型单一固定端元也行。3.3 结果写出与投影转换GeoTIFF 保存和 Albers 重投影算完 FVC 要写出 GeoTIFF保持和原数据一致的投影和范围。如果要做面积统计再重投影到 Albers。# 写出 FVC 栅格保持原投影 profile.update(dtypefloat32, nodatanp.nan, count1) with rasterio.open(FVC_2015_China_1km.tif, w, **profile) as dst: dst.write(fvc.astype(float32), 1) # 重投影到 Albers 等积投影命令行方式更稳 # gdalwarp -t_srs projaea lat_125 lat_247 lat_00 lon_0105 datumWGS84 unitsm \ # FVC_2015_China_1km.tif FVC_2015_China_Albers.tif逻辑说明profile.update把 dtype 改成 float32nodata 改成 nancount 保持 1。写出后用gdalwarp重投影Albers 参数里lat_1和lat_2是中国常用双标准纬线lon_0105是中央经线。重投影后像元大小会变面积统计更准。参数说明gdalwarp的-t_srs后面跟目标投影字符串-r可以指定重采样方法默认 near 适合分类数据连续数据用 bilinear 或 cubic。FVC 是连续值建议-r bilinear。4. 避坑与排查2015 年 NDVI 数据处理的 5 个血泪教训4.1 现象全国均值算出来只有 0.2明显偏低原因NoData 值没掩掉比如 -3000 被当成有效值参与平均。解决先跑gdalinfo确认 NoData再用ndvi ! nodata掩膜。如果 NoData 是 nan用np.isnan判断。4.2 现象FVC 出现大量负值或大于 1原因NDVI_soil 和 NDVI_veg 取反了或者分位数算在了包含 NoData 的数组上。解决确保分位数计算前已经掩掉无效像元并且 NDVI_veg 一定大于 NDVI_soil。如果还是溢出用np.clip截断。4.3 现象重投影后面积统计和原始差很多原因正弦投影下高纬度像元面积被拉伸直接统计会偏大。解决先重投影到 Albers 等积投影再统计或者用rasterio的xy方法逐像元算实际面积。4.4 现象QA 波段没读云污染没剔除原因只用了 NDVI 波段忽略了 QA 里的云标志。解决如果数据带 QA 波段按位解析把云、雪、阴影对应的像元掩掉。常见做法是qa 0b11 ! 0判断云。4.5 现象内存爆了程序被 kill原因全国 1km 栅格虽然不大但中间变量太多float64 转 float32 没做。解决读数据时直接astype(float32)中间计算用np.nanmean而不是先填充再算及时del不用的变量。5. 进阶技巧用 2015 年 NDVI 做年际对比与趋势验证如果你手上不止 2015 年一年而是有 2010-2020 的序列那这份 2015 年数据就是基准年。我一般会做两件事一是算 2015 年相对 2010 年的 NDVI 变化率二是用 Theil-Sen 斜率做趋势检验。Theil-Sen 比线性回归稳对异常值不敏感适合遥感数据。from scipy.stats import theilslopes # 假设 ndvi_stack 是多年 NDVI 数组形状 (年份, 行, 列) # 逐像元算 Theil-Sen 斜率 slope np.full(ndvi_stack.shape[1:], np.nan, dtypefloat32) for i in range(ndvi_stack.shape[1]): for j in range(ndvi_stack.shape[2]): ts ndvi_stack[:, i, j] if np.isnan(ts).any(): continue slope[i, j] theilslopes(ts, np.arange(len(ts)))[0]这段代码慢但胜在稳。如果嫌慢可以用xarray的polyfit或者并行化。验证方法把 2015 年单年 FVC 和多年趋势叠加看哪些区域是「高覆盖但下降」或「低覆盖但上升」这些往往是生态工程的重点区。注意Theil-Sen 要求时间序列至少 5 年以上少于 5 年噪声太大别硬跑。我自己的习惯是拿到任何一年的 NDVI 数据先跑一遍gdalinfo再算有效像元占比然后才动 FVC。2015 年这份数据我前后用过三次第一次没掩 NoData全国均值 0.18差点把结论写反第二次忘了重投影面积统计多了 7%第三次才顺过来。希望帮到你。本文还有配套的精品资源点击获取
MODIS 2015年中国1km NDVI数据：从读取到FVC反演的Python实战
简介本资源为2015年中国区域1km分辨率植被指数NDVI空间分布数据集面向遥感、生态、农业及地理信息领域的科研人员与学生可用于植被覆盖监测、土地退化评估、气候变化分析等场景。数据源自NASA MODIS系列MOD13A3产品经子数据集提取、拼接、投影转换、单位换算与裁剪后采用最大合成法生成年度NDVI投影为Albers等面积圆锥投影椭球WGS84中央经线105°标准纬线25°与47°。压缩包共5个文件约20.01MB包含1个tif栅格主文件、1个tfw坐标文件、2个xml元数据及1个txt说明文档便于直接加载与坐标配准。已有312人学习下载。读者可获取完整年度NDVI栅格用于空间分布制图、时序对比与模型输入并借助说明文档快速理解数据来源与处理流程。1. 从一份 2015 年中国 1km NDVI 数据说起它到底能干什么如果你手上正好有一份 MODIS 2015 年中国 1km 植被指数NDVI空间分布数据集第一反应大概率不是「这数据真漂亮」而是「我拿它能算什么」。我最早接触这类数据是为了做县域尺度的植被覆盖度FVC反演当时手头只有行政边界和一堆气象站点最缺的就是一张能覆盖全国、时间对得上、分辨率又不至于太粗的底图。NDVI 正好补上了这个缺口——它用近红外和红光波段的反射率差值归一化直接反映植被绿度1km 像元意味着全国大约 960 万平方公里被切成不到一千万个格子单省跑一遍内存扛得住做趋势分析也不至于被云污染撕得七零八落。这份数据集的核心价值在于「空间分布」四个字。它不是单点观测不是站点插值而是栅格化的面状覆盖能直接和土地利用、气候区划、DEM 做叠加。适合谁用做生态评估的、算 FVC 的、盯退耕还林效果的、写区域规划报告的甚至只是想给论文配一张全国植被长势图的都能从这份 2015 年的切片里挖出东西。但前提是你得先搞清楚它的投影、填充值、有效范围否则后面每一步都是坑。2. MODIS NDVI 的物理含义与 2015 年 1km 产品的选型逻辑2.1 NDVI 为什么能代表植被从波段反射到归一化差值NDVI 的公式简单到一行就能写完(NIR - Red) / (NIR Red)。但简单不代表随便用。健康植被的叶肉组织强烈反射近红外同时叶绿素吸收红光所以 NIR 高、Red 低差值大NDVI 就高。裸土、水体、云、雪则相反NDVI 趋近于零甚至为负。MODIS 传感器在 2015 年已经稳定运行了十几年波段响应函数成熟辐射定标精度在 2% 以内这是它比很多国产传感器更早被拿来做长时序分析的原因。但要注意NDVI 不是万能的。它饱和得早——当叶面积指数超过 4 左右NDVI 就涨不动了所以密林区容易低估变化。另外它受土壤背景影响大稀疏植被区信号弱。2015 年这个时间点选得也巧既避开了 2000 年代初传感器刚上天时的数据波动又没到后来部分波段衰减明显的阶段做基准年比较合适。2.2 为什么是 1km 而不是 250m 或 500m分辨率与全国覆盖的权衡MODIS 的 NDVI 产品有 250m、500m、1km 三档。250m 最细但全国跑一遍数据量直接翻四倍而且 250m 只有红光和近红外两个波段云掩膜和大气校正的辅助信息少质量控制字段不如 1km 丰富。500m 是折中但很多做省级以下单元分析的人会发现 500m 在县级尺度上还是太粗一个县没几个像元统计意义弱。1km 产品的优势在于第一它带了详细的 QA 波段能逐像元判断云、雪、阴影、气溶胶影响第二全国范围单年数据量在几百 MB 到 1GB 量级普通笔记本能处理第三很多已有的土地利用、气象、土壤数据集本身就是 1km 或更粗匹配起来不用重采样。我一般会建议做全国或大区域趋势分析1km 足够做地块级或城市内部绿化评估直接上 250m 或更高分辨率别在 1km 上硬抠。2.3 2015 年中国区数据的获取与格式确认这份数据集常见格式是 GeoTIFF 或 HDF投影多为正弦投影Sinusoidal也有已经转成 Albers 等积投影的版本。拿到手第一件事不是打开看而是用gdalinfo把元数据读出来。下面这段命令是我每次拿到新栅格必跑的gdalinfo MODIS_NDVI_2015_China_1km.tif重点看几个字段Size is后面的行列数Coordinate System is里的投影和基准面NoData Value填充值Band 1的Type和STATISTICS。如果 NoData 是 -3000 或 32767 这类值后面计算前必须掩掉否则全国均值会被拉低一大截。投影如果是正弦的做面积统计前先重投影到等积投影不然高纬度地区面积会失真。提示不要直接用 QGIS 打开就目测。先跑gdalinfo把 NoData、投影、行列数记在脚本注释里后面每一步都引用这些参数。3. 用 Python 把 2015 年 NDVI 栅格跑通从读取到 FVC 反演3.1 读取栅格与有效像元掩膜rasterio 最小代码读取用rasterio最稳它比 GDAL 的 Python 绑定更友好。下面这段代码做了三件事读数据、读 NoData、生成有效掩膜。import rasterio import numpy as np path MODIS_NDVI_2015_China_1km.tif with rasterio.open(path) as src: ndvi src.read(1).astype(float32) # 读第一波段 nodata src.nodata # 读取 NoData 值 profile src.profile # 保存元数据供后续写出 print(行列数:, src.height, src.width) print(投影:, src.crs) print(NoData:, nodata) # 构建有效掩膜排除 NoData 和 NDVI 异常值 valid np.ones(ndvi.shape, dtypebool) if nodata is not None: valid (ndvi ! nodata) valid (ndvi -1.0) (ndvi 1.0) # NDVI 理论范围 ndvi_valid np.where(valid, ndvi, np.nan) print(有效像元占比: {:.2f}%.format(valid.mean() * 100))逻辑说明src.read(1)读的是第一波段如果文件里还有 QA 波段波段顺序要确认。nodata可能是 None所以先判断。NDVI 理论范围是 -1 到 1超出这个范围的像元基本是填充或异常直接掩掉。最后用np.nan替换无效值方便后面用np.nanmean统计。参数说明astype(float32)是为了省内存全国 1km 栅格大概 8000×6000 左右float32 占不到 200MB。如果内存紧张可以分块读但 2015 年单年数据一般不用。3.2 从 NDVI 到植被覆盖度 FVC像元二分模型参数怎么定FVC 的像元二分模型公式是FVC (NDVI - NDVI_soil) / (NDVI_veg - NDVI_soil)。关键在于 NDVI_soil 和 NDVI_veg 怎么取。常见做法是取累计频率的 5% 和 95% 分位数而不是固定值。下面代码演示# 基于有效像元计算分位数 ndvi_flat ndvi_valid[~np.isnan(ndvi_valid)] ndvi_soil np.percentile(ndvi_flat, 5) # 裸土端元 ndvi_veg np.percentile(ndvi_flat, 95) # 纯植被端元 print(NDVI_soil {:.4f}, NDVI_veg {:.4f}.format(ndvi_soil, ndvi_veg)) # 像元二分模型计算 FVC fvc (ndvi_valid - ndvi_soil) / (ndvi_veg - ndvi_soil) fvc np.clip(fvc, 0, 1) # 截断到 0-1 # 统计全国均值 print(2015 年全国 FVC 均值: {:.4f}.format(np.nanmean(fvc)))逻辑说明分位数取 5% 和 95% 是经验值目的是排除极端异常像元。如果研究区是纯森林NDVI_veg 会偏高如果是干旱区NDVI_soil 可能接近 0。所以最好分区域定参数不要全国一刀切。np.clip把超出 0-1 的值截断避免负 FVC 或大于 1 的荒谬值。参数说明分位数可以调比如 2% 和 98% 更激进5% 和 95% 更保守。我一般先跑一遍看直方图如果 NDVI 分布双峰明显分位数法很稳如果单峰说明研究区植被类型单一固定端元也行。3.3 结果写出与投影转换GeoTIFF 保存和 Albers 重投影算完 FVC 要写出 GeoTIFF保持和原数据一致的投影和范围。如果要做面积统计再重投影到 Albers。# 写出 FVC 栅格保持原投影 profile.update(dtypefloat32, nodatanp.nan, count1) with rasterio.open(FVC_2015_China_1km.tif, w, **profile) as dst: dst.write(fvc.astype(float32), 1) # 重投影到 Albers 等积投影命令行方式更稳 # gdalwarp -t_srs projaea lat_125 lat_247 lat_00 lon_0105 datumWGS84 unitsm \ # FVC_2015_China_1km.tif FVC_2015_China_Albers.tif逻辑说明profile.update把 dtype 改成 float32nodata 改成 nancount 保持 1。写出后用gdalwarp重投影Albers 参数里lat_1和lat_2是中国常用双标准纬线lon_0105是中央经线。重投影后像元大小会变面积统计更准。参数说明gdalwarp的-t_srs后面跟目标投影字符串-r可以指定重采样方法默认 near 适合分类数据连续数据用 bilinear 或 cubic。FVC 是连续值建议-r bilinear。4. 避坑与排查2015 年 NDVI 数据处理的 5 个血泪教训4.1 现象全国均值算出来只有 0.2明显偏低原因NoData 值没掩掉比如 -3000 被当成有效值参与平均。解决先跑gdalinfo确认 NoData再用ndvi ! nodata掩膜。如果 NoData 是 nan用np.isnan判断。4.2 现象FVC 出现大量负值或大于 1原因NDVI_soil 和 NDVI_veg 取反了或者分位数算在了包含 NoData 的数组上。解决确保分位数计算前已经掩掉无效像元并且 NDVI_veg 一定大于 NDVI_soil。如果还是溢出用np.clip截断。4.3 现象重投影后面积统计和原始差很多原因正弦投影下高纬度像元面积被拉伸直接统计会偏大。解决先重投影到 Albers 等积投影再统计或者用rasterio的xy方法逐像元算实际面积。4.4 现象QA 波段没读云污染没剔除原因只用了 NDVI 波段忽略了 QA 里的云标志。解决如果数据带 QA 波段按位解析把云、雪、阴影对应的像元掩掉。常见做法是qa 0b11 ! 0判断云。4.5 现象内存爆了程序被 kill原因全国 1km 栅格虽然不大但中间变量太多float64 转 float32 没做。解决读数据时直接astype(float32)中间计算用np.nanmean而不是先填充再算及时del不用的变量。5. 进阶技巧用 2015 年 NDVI 做年际对比与趋势验证如果你手上不止 2015 年一年而是有 2010-2020 的序列那这份 2015 年数据就是基准年。我一般会做两件事一是算 2015 年相对 2010 年的 NDVI 变化率二是用 Theil-Sen 斜率做趋势检验。Theil-Sen 比线性回归稳对异常值不敏感适合遥感数据。from scipy.stats import theilslopes # 假设 ndvi_stack 是多年 NDVI 数组形状 (年份, 行, 列) # 逐像元算 Theil-Sen 斜率 slope np.full(ndvi_stack.shape[1:], np.nan, dtypefloat32) for i in range(ndvi_stack.shape[1]): for j in range(ndvi_stack.shape[2]): ts ndvi_stack[:, i, j] if np.isnan(ts).any(): continue slope[i, j] theilslopes(ts, np.arange(len(ts)))[0]这段代码慢但胜在稳。如果嫌慢可以用xarray的polyfit或者并行化。验证方法把 2015 年单年 FVC 和多年趋势叠加看哪些区域是「高覆盖但下降」或「低覆盖但上升」这些往往是生态工程的重点区。注意Theil-Sen 要求时间序列至少 5 年以上少于 5 年噪声太大别硬跑。我自己的习惯是拿到任何一年的 NDVI 数据先跑一遍gdalinfo再算有效像元占比然后才动 FVC。2015 年这份数据我前后用过三次第一次没掩 NoData全国均值 0.18差点把结论写反第二次忘了重投影面积统计多了 7%第三次才顺过来。希望帮到你。本文还有配套的精品资源点击获取
