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 env
sourceDir = 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 + '\\' + eviName
print(outname) #打印文件名用于检查
arcpy.ExtractSubDataset_management(hdf, outname)#可以改为(hdf, outname,0),0表示第一层,即NDVI,具体各层是啥可以通过arcgis的ExtractSubDataset工具查看。不设置就默认提取第0层
使用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%左右,
参考他人论文研究成果,确定那种结果最可靠。