使用MOD13A1产品NDVI数据计算植被覆盖度

MOD13A1数据分辨率500m，相比A2与A3的分辨率略高，计算结果可能更准确。也要看研究区大小，我的研究区较大，且其他数据分辨率是1KM，其实没必要使用500m的分辨率。
MOD13Q1分辨率250m，适合更小的研究区使用。

计算过程：

1、下载数据：

https://search.earthdata.nasa.gov/
NASA官网下载，推荐使用批量下载工具：DOWNthemALL
我这里下载的是2020年的NDVI数据

2、使用python或ARCGIS将HDF转TIF

arcgis可以直接打开HDF数据，但是比较慢，使用python对HDF数据进行格式转换：

# coding=utf-8
import os
import arcpy
from arcpy import envsourceDir = r'E:\GISdata\MOD13A1\2020HDF'  # 输入
targetDir = r'D:\MOD13A1TIF'  # 输出
#检查Spatial —ArcGIS Spatial Analyst 扩展模块是否许可
arcpy.CheckOutExtension("Spatial")
#设置工作空间
env.workspace = sourceDir
arcpy.env.scratchWorkspace = sourceDir
#读取栅格列表
hdfList = arcpy.ListRasters('*', 'hdf')
# for 循环
for hdf in hdfList:# print hdf  #打印文件名用于检查# 替换文件名后缀eviName = os.path.basename(hdf).replace('hdf', 'tif')outname = targetDir + '\\' + eviNameprint(outname)   #打印文件名用于检查arcpy.ExtractSubDataset_management(hdf, outname)

使用arcpy转hdf到tif的好处是，直接生成了投影文件等，再在arcgis中打开时，不需要建立影像金字塔等。

3、使用最大值合成法进行镶嵌

加载到arcgis中发现影像的范围为-2000到9xxx。

官网显示数据的有效区间是（-2000——10000）,扩大了10000倍

建模进行镶嵌，选择最大值合成，如果电脑处理速度还行，可以直接改成32bit浮点数，要是速度不行处理完再改浮点数类型

需要处理一段时间

按说这样处理的逻辑没问题，可是出来的图像是不完整的。
———————————————————————————————————————
使用另一种方法，直接使用镶嵌至新栅格工具：

在导入的时候只需要按住shift，全选就可以一次把几十个都导入，不需要一个个选择。这样处理的结果是争取的。

4、转换单位并掩膜提取

使用掩膜提取工具进行掩膜提取，之后使用栅格计算器，对影像除以10000。我不是很明确是否要把小于0的数据直接让其等于0，所以在后面的操作中产生了两个版本，一个是没有将小于0的进行转换，另一个是将小于0的值都让其等于0（remove negative value）。

通过将两幅图像相减，发现我的研究区内，NDVI是负值的栅格其实很少很少，只有很个别的栅格，所以我认为很多情况下，是否将负值转为0的影响不大。

5、计算植被覆盖度：

有的人认为，NDVIsoil与NDVIveg在NDVI值的范围是0-1的情况下，可以直接使NDVIsoil=NDVImin, MDVIveg =NDVImax。但不适用我现在的情况。

对于置信度的选择，有的研究选择的5%，这里的研究选用的0.5%，在此选用5%与1%置信度的方法进行对比，看结果的差异。

不太懂如何从arcgis中查看置信度，所以使用累计百分比代替置信度。
打开影像的符号系统，选择使用分位数进行分类。

下图为累计百分比为5%的去掉负值的情况：

去掉负值的NDVI的累计百分比5%：0.229845，95%:0.896008
没有去掉负值的NDVI累计百分比5%：0.228506，95%：0.893235
去掉负值的NDVI的累计百分比1%：0.124662，99%:0.942756
没有去掉负值的NDVI累计百分比0.5%：0.128341，99%：0.943318

分别令累计百分比1%或5%作为NDVIsoil，95%或99%为NDVIveg，计算植被覆盖度。

1%累计百分比的计算结果比5%的植被覆盖度大2%左右，

参考他人论文研究成果，确定那种结果最可靠。