前言：

2023.03.07修改：HEG批量拼接处理代码缩进问题没改完……

2022.12.30修改：补Python编程基础的过程中发现使用Print也可以向文件里写东西，查看写过的代码，然后发现HEG批量拼接处理代码中有的地方在复制粘贴至博客的时候出现了缩进错误,评论区沉冤昭雪……，已修改。

#关于print写文件
如：
f = open('test.txt','a+')
print('Hello World!',file=f)
f.close()

2021.12.25修改：给目录添加了一个下载数据，纠正了python栅格计算误区，其实把变量定义为浮点数就可以解决溢出的问题（我是笨蛋）。

虽然百度上啥都有，但是现在垃圾广告和复制粘贴的垃圾网站太多，有用的信息基本都被淹没了，前后历时三四天总算弄明白了MODIS数据怎么下载，HEG数据的批量处理，使用python进行栅格计算的一些事情，就在这儿写一写。

使我感触最深的是：作为一名程序员，首先你要学会百度，不是百度全名，而是想好步骤按步骤百度（因为具体问题要具体分析）

我弄明白的只是一种方法。

数据来源：Find Data - LAADS DAAC

NASA官网的MOD13A1数据集，500M空间分辨率

数据选择阶段需要VPN，下载不需要

目录

数据下载： 

HEG安装：

批量使用HEG：

批量掩膜提取：

栅格值计算：

---

数据下载： 

---

我推荐你使用Google浏览器，因为它可以全文翻译页面，而且翻译效果也不是很差。

点红圈里的按钮就可以

如果你发现你的浏览器那里没有那个按钮也不要紧，在页面空白处右键单击一下，你会看到右键菜单里有“翻成中文（简体）”这个选项，很简单对吧。

然后是下载数据的事情，右侧有6个大按钮，后5个在你下订单（免费，不花钱）之前是没有用的，你不用管，单击“按产品搜索”就会进入我截图的这个页面。

你应该也看到页面的上方有5个步骤，分别是

产品（products） 时间（time） 地点（location） 档案（files） 查看订单（review & order）

那咱就一步一步来呗

如果你是英文状态还看不懂英文的话，你可能要问产品说明在哪儿，我怎么找到我需要的东西，但是中文就不需要烦恼，因为上面写着

那咱们就随便选一个进行下一步，

值得一提的是，他们考虑到了你一次想下载多种数据的情况

然后是日期选择，也是个挺折磨人的地方，不知道为什么，日历里面选择的日期并不能出现在文字框里面，也许是我的问题。

个人推荐按文本框里面的格式直接输入日期，然后点下方的添加日期

与产品相同，这个你也可以选择好几个时间段，同时右边的X可以删除你选错的时间段，时间段选好之后就进入了地点选择

如果你要下个全球或者某个国家的那再好不过了，省去好多功夫。

网站允许你以国家，瓷砖（地图分的块，我不知道为什么Google把它翻译为瓷砖），自己画个框，和输入经纬度的方式选择区域。

国家和瓷砖我就不说了，这里会显示你鼠标的经纬度坐标，所以你至少需要先找到你研究区域所跨的经纬度，百度或者先下载一幅世界或者国家的影像在arcmap里面瞅瞅都行，然后拉个框

输入坐标也是一种方法，但是你要注意，不需要输入WNES，也不需要输入°，

坐标系统采用的是经度-180到180，维度90到-90，按西经，北纬，东经，南纬的顺序依次输入，数值之间用英文逗号分隔，看我输入的样例。

之后就是选择文件了，当然是全选了，你以为呢？

查看订单的时候，你需要一个账户，没有的话创建一个，也就是点登记（register），注册并不困难，只需要一个邮箱，建议不要用QQ邮箱，会出问题。

注册登记后就回到了订单页面，你可以继续添加其他数据（add another search），或者下订单（submit order），下订单之后不要着急，等着官方给你发订单已处理邮件，大概两三分钟。

收到邮件后等个两三分钟左右，因为网站数据更新会延迟，点击右侧的过往订单（past orders）

订单会为你保存一周左右，你要是不下载它就帮你删除了，就那些白色的

然后可用订单是绿色的，就是我刚刚弄得实例文件，红圈里面有个订单号，点击后会跳转下载页面。

啊，下图就是没有等待出现的错误，你的数据还没有载入到你的用户里面，网站找不到它们，可以刷新，也可以后退，都行

全选，然后下载，当然，你可以先下载一幅看看行不行

然后是个很重要的事情

单次下载文件不要超过二百个，不然它会崩溃，想体验可以试一试，当连续下载文件数接近200的时候，网站要么卡死不下载，要么给你下载重复的文件或者跳着给你下载，给下载后的数据校对带来极度不适。

数据下载完之后会有一个以checksums开头的文件，帮助你校对自己的文件少了没有。

下载过程建议打开VPN，虽然不使用VPN的流量，但是关闭VPN有时会出现下载文件错误的情况。

至此我们完成了数据下载，hdf格式的文件。

HEG安装：

人家的数据肯定要用人家的处理工具，（我十分佩服水经注，80GB的影像文件竟然可以压缩成2GB的二进制文件，实在厉害）

显然HEG并不是一个冷门的东西，因为我在打不开其他作者给的下载链接的时候，直接Google搜索了HEG，然后就得到了如下页面（翻译成中文）

贴个链接：

HEG: HDF-EOS to GeoTIFF Conversion Tool - Data Access Services - Earthdata Wiki

点下载，然后选一个

从上至下依次为，HEG用户指南（非常详细的PDF，如果你的英文够好，看它吧，别再看这个文章了），2.15注意事项，Linux版HEG，MAC版HEG，Windows版HEG，HEG安装指南。

我想你会喜欢这个页面的，也是从Google得到的，

我觉得世界上的网站排布好迷啊，我总是找不到我需要的功能块在哪儿，也许是我的问题。

这个页面算是个简单版带图片的用户指南，，划红线的链接点击也可以下载HEG用户指南

贴个页面链接：HDF-EOS to GeoTIFF Conversion Tool (HEG)

之后就是安装的事儿了

软件需要JAVA运行环境的支持，所以没有安装它们需要先装一个（就像我）

Java很容易找到，百度JDK和Java都可以

就在官网下载一个，然后安装，HEG好像需要Java1.8及以上版本的支持，如果你的版本低，记得更新（因该不会出现这情况，一般电脑上软件版本越老说明这个人越厉害，只有我这种萌新才会去装最新的，因为我不知道哪个好用，于是导致了之后的一些悲剧）

安装的时候注意一下你的安装路径，之后需要用到。

贴个java下载链接

https://download.oracle.com/java/17/latest/jdk-17_windows-x64_bin.exe

然后是安装HEG，解压会得到一堆文件，建议你建个文件夹再解压，双击install.bat就可以安装，

它会询问你要不要安装，简直是废话

给它定义安装路径，我C盘装了一个，这次装D盘，顺带提一句，如果你有机械硬盘，千万别把任何软件装机械硬盘里面，固态和机械的存取速度差几百倍，可不是一般的影响性能。

大多数国外软件对中文路径很敏感，中文路径经常会弄死它们，你可以给文件夹起个拼音名字，然后继续y。

这里就是需要你记路径的原因了，找到你的Java文件夹里的bin，写入路径

我这里犯了个错误，不知道为什么，java和python语言里面，路径的一般用的是斜杠‘/’而不是反斜杠‘\’，而Windows默认路径名字就是反斜杠，需要改一下，不然会错

如果你的路径中有Program files，你需要注意，两个单词之间有个空格，别打空格或者空格前加一个转义字符‘\’也行。

虽然错误不会影响你的安装，但是会影响你的使用，

接下来输入用户名和密码，其实没什么用，就安装成功了

安装成功后你会在\HEG\HEG_Win\bin 目录下找到一个文件

HEGTool.bat

双击就可以看到GUI，可视化的操作窗口，第一列是供你查看你导入的数据，选择你需要的那一部分，第二列选择输出路径，重采样方法，投影坐标系，第三列可以保存（Save）参数文件或者执行（Run），操作很简单，我就不说那么多了

 

批量使用HEG：

写给不看和看不懂说明书的盆友（比如我）

如果你只是转tif，不需要拼接，那太好了，打开选择第二个，找到数据文件夹，选择你的第一个hdf，选好你需要的东西和输出路径，别点Run，点Batch Run，其他交给它，傻瓜式操作。

使用python脚本批量调用主要适用于拼接的情况，

首先点击tool→stitch/subset转换到拼接工具，拼接和格式转换以及坐标变化可同时进行，你并不需要搞两遍，很Nice不是吗。

在开始之前你需要知道一些东西，比如MODIS数据的命名规则：

数据分为五块：MOD13A1成品来源，A2013113是影像数据日期编码即2013年第113天，h25v04是地图分幅号，006.……是处理数据的日期编码，hdf文件格式。

也就是说，你需要把同一天的拼一起

先导入数据，点第一个，一个一个导入，不要点open下拉菜单的第二个，我曾天真的以为它功能很完善，后来发现我错了，也可能是我太笨没学会怎么用。我就不全导入了。

我选了NDVI，browse选了自己指定的文件夹，其他默认，Accept，之后不要点运行，点Save！

保存后的是一个txt文件，Unix编码（这一点很重要）的txt文件，可以用记事本打开

这里记录了处理参数，然后咱们需要做的就是批量更改参数并调用HEG来处理。

需要修改的有这几个参量

NUMBER_INPUTFILES 文件个数

INPUT_FILENAMES 文件名（多个文件名之间以‘|’分隔）

OUTPUT_FILENAME 输出文件名

OUTPUT_STITCHED_FILENAME 这一行属于过程文件的输出路径，虽然选了NO，但是批量调用的时候它还是会给你输出，建议这一行删掉（不影响程序运行）

会使用python编程的可以走了

然后我来讲讲我作为一个完全没接触过python的人写出来的东西，我建议使用notepad++编写python文件，很好用。Python对空格很敏感，一定要注意。

路径要使用‘/’连接

# -*- coding: utf-8 -*-

"""

调用HEG相关工具批量完成MODIS数据拼接工作。

"""

import os


# 设置HEG相关环境变量，简单来说就是把data，TOOLKIT_MTD,bin的路径告诉系统，方便cmd执行，别问我为什么，用户手册上写的

os.environ['MRTDATADIR']='C:/Works/HEG/HEG_Win/data'

os.environ['PGSHOME']='C:/Works/HEG/HEG_Win/TOOLKIT_MTD'

os.environ['MRTBINDIR']='C:/Works/HEG/HEG_Win/bin'



#路径，程序，输出，输入都得写，一个for循环获取一下需要被处理的所有文件

# 设置HEG的bin路径

hegpath = 'C:/Works/HEG/HEG_Win/bin'

# 指定处理模块的可执行程序文件subset_stitch_grid，

hegdo = os.path.join(hegpath, 'subset_stitch_grid.exe')

hegdo = hegdo.replace('\\', '/') # 全路径以“/”连接


# 指定输入数据的路径

inpath = r'C:\Users\mrlon\Desktop\share\modis\2021'

inpath = inpath.replace('\\', '/')

# 指定输出数据的路径

outpath = r'C:\Users\mrlon\Desktop\share\modis\2021_out'

outpath = outpath.replace('\\', '/')

# os.chdir(inpath) #改变当前工作目录到输入数据目录


# 获取当前文件夹下的所有hdf文件

allfiles = os.listdir(inpath)

allhdffiles = []

for eachfile in allfiles:

    if os.path.splitext(eachfile)[1] =='.hdf': #这一句是为了把不是hdf的文件给排除，可以去掉

        allhdffiles.append(eachfile)

#用print把文件名显示一下，也可以去掉

print('--'*20)

print('文件数量为:', len(allhdffiles),',所有hdf文件如下')

print('  '+'\n  '.join(allhdffiles))

print('--'*20)



# 需要在HEG工具中生成一个参考的文件，咱已经生成了一个，设置文件存储路径

prmpath = r"C:\Users\mrlon\Desktop\123 \ 1_gridstitch"

prmpath = prmpath.replace('\\', '/')

#然后就是一个循环，需要多少个参数文件，咱就循环多少次写多少个
'''
2022.12.30 我的代码应该是一个单次循环完成一个txt文件
所以prmfilename=prmpath +'/'+ allhdffiles[i] +'_gridstitch'
以及fo.close()
是循环内的语句，然后复制粘贴过程少了缩进
'''

i=0

while i < len(allhdffiles):

    #变量prm暂存过会儿需要写进txt的文本

    prm=[

        'NUM_RUNS = 1\n',

         'BEGIN\n',

         'NUMBER_INPUTFILES = 7\n', #提前根据命名规则看一下你需要把多少个文件合在一起，这里就写几

        'INPUT_FILENAMES='+inpath+'/'+allhdffiles[i]+'|'+inpath+'/'+allhdffiles[i+1]+'|'+inpath+'/'+allhdffiles[i+2]+'|'+inpath+'/'+allhdffiles[i+3]+'|'+inpath+'/'+allhdffiles[i+4]+'|'+inpath+'/'+allhdffiles[i+5]+'|'+inpath+'/'+allhdffiles[i+6]+'\n', #这里是路径，我的i从0开始，一次循环加7，于是我这里就有七个文件名；

         'OBJECT_NAME = MODIS_Grid_16DAY_500m_VI|\n',

         'FIELD_NAME = 500m 16 days NDVI|\n',

         'BAND_NUMBER = 1\n',

         'SPATIAL_SUBSET_UL_CORNER = ( 49.999999996 78.324437349 )\n',

         'SPATIAL_SUBSET_LR_CORNER = ( 29.999999997 155.572382658 )\n',

         'OUTPUT_OBJECT_NAME = MODIS_Grid_16DAY_500m_VI|\n',

         'OUTGRID_X_PIXELSIZE = 0.0041666733029590754\n',

         'OUTGRID_Y_PIXELSIZE = 0.0041666733029590754\n',

         'RESAMPLING_TYPE = BI\n',

         'OUTPUT_PROJECTION_TYPE = GEO\n',

         'ELLIPSOID_CODE = WGS84\n',

         'OUTPUT_PROJECTION_PARAMETERS = ( 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0  )\n',

         'OUTPUT_FILENAME = '+outpath+'/'+allhdffiles[i]+'_out.tif\n',

#输出文件名，路径+名字+后缀名，你可以随意命名

         'SAVE_STITCHED_FILE = NO\n',

         #'OUTPUT_STITCHED_FILENAME = '+outpath+'/'+allhdffiles[i]+'_stitched.hdf\n',

#你应该看到我给这里多加了个#号，对，我把它注释掉了，不然他会给我生成一堆临时文件，而且还挺大

         'OUTPUT_TYPE = GEO\n',

         'END\n',

        

         ]

    prmfilename=prmpath +'/'+ allhdffiles[i] +'_gridstitch'

#设置参数文件输出路径和名字

    i+=7

    prmfilename=prmfilename.replace('\\', '/')

    #这里一定要注意，设定换行符为‘\n’,否则由于在windows系统下默认换行符为‘\r\n’,则无法运行成功

    fo=open(prmfilename,'w',newline='\n')

    fo.writelines(prm)

    fo.close()

#写入参数文件

   

#第一个循环结束



#第二个循环，批量调用并执行，七个文件，j+7

j=0

while j < len(allhdffiles):

    prmfilepath=prmpath +'\\'+ allhdffiles[j] + '_gridstitch'

#参数文件路径

#2023.03.07 修改了上方语句的缩进问题

    prmfilepath=prmfilepath.replace('\\', '/')

    try:

        resamplefiles = '{0} -P {1}'.format(hegdo, prmfilepath)

        os.system(resamplefiles) 

#这两句就是调用cmd，固定写法，hegdo、prmfilepath是我定义的变量，分别储存程序和参数文件路径     

        print(allhdffiles[j] + ' has finished')

        j+=7

    except:

        # 提示错误信息

        print(eachhdf + 'was wrong')

 

      

完事儿，哦不对，还有运行的事儿

我再次推荐使用notepad++，因为HEG的参数文件是Unix编码的txt，而用python直接执行的话是Windows的txt，HEG识别不了，会出错，notepad++可以一键转换，转换后run就可以

初次运行需要你输入一下你的python.exe在哪儿

批量掩膜提取：

对完全不懂python的我，这应该是最简单的部分了，esri永远的神！

首先你需要准备研究区边界的shpfile；

我觉得没有arcmap的人应该也不会看这篇文章，所以就不写安装过程了，需要的可以百度。

打开模型构建器：

插入，迭代器，栅格：

选择tif存储的文件夹目录，递归适用于一个文件夹里还有数个子文件夹的情况。

在工具箱找到你需要的工具，这里我需要掩膜，就把掩膜拖进去，

检查一下拓展模块，spatial analyst需要勾选，不然你用不了

拖入掩膜，用链接工具连接一下：

在输出文件上点右键，打开，找好输出文件夹，输出名称设置为“%名称%”

这样可以保证你输出的结果名称和裁剪前名称一样，所以最好不要放在同一个文件夹。

另外注意输出路径不要有中文，不然会报错。

然后run就可以。

当你跑第二遍却发现只提取了一个文件的时候，打开迭代器，随便点点，比如递归那里打个勾再取消，然后确定，运行就好，因为迭代器偶尔会出现到最后一个文件不能自动重置的情况。

栅格值计算：

这里就又回到了python，虽然用arcmap打开可以直接获取均值，极值等东西，但是一个个打开也太难了吧，我有九千多个文件啊

于是使用python脚本进行栅格计算就好得多了

正式开始前你需要打开一张你的tif，在它的属性→源，查看像素数量，均值，Nodata值，极值，便于你去除无用的像素点和后期校对数据

另外你需要注意，同源数据的Nodata值基本一致，但偶尔会出现不一致的情况，当你发现结果很奇怪的时候，记得回到arcmap去看一眼。

然后就是编码：

#这个程序获取TIF文件的所有像元值做平均并将结果写入Excel表格。也就是打开栅格文件，遍历所有像元，像元值做和，求商，保存

import os #导入os包

from osgeo import gdal #导入gdal包，我忘了什么用，好像不需要

#导入numpy包（支持高维数组和矩阵运算，也提供了许多数组和矩阵运算的函数）

import numpy as np

#导入openpyxl包用来搞excel

from openpyxl import Workbook

import openpyxl



#指定工作文件夹，虽然我想用三层循环，也就是文件夹里套文件夹再套文件夹一运行全部计算，后来发现数据量太大这样做有很大风险，于是就套了两层，首先获取陕西文件夹里面的所有文件夹

#路径

inpath_province = r'C:\Users\mrlon\Desktop\share\city_extract\陕西'

inpath_province = inpath_province.replace('\\','/')

#获取全部文件夹

allfiles_province = os.listdir(inpath_province)

#定义一个列表，装列表里，数组也行

workfiles_province = []

for eachfile in allfiles_province:

    workfiles_province.append(eachfile)

#传统的print，看看对不对

print('--'*20)

print('文件数量为:', len(workfiles_province),',所有文件夹如下')

print('  '+'\n  '.join(workfiles_province))

print('--'*20)



#主体，一个for循环，w代表的是文件夹的个数

w=0

while w<len(workfiles_province):

 n=0

 # 指定输入tif数据的路径，w=0先从第一个文件夹开始呗

 inpath = inpath_province +'\\'+workfiles_province[w]

 inpath = inpath.replace('\\', '/')

# 获取当前文件夹下的所有tif文件写到列表里

 allfiles = os.listdir(inpath)

 allhdffiles = []

 for eachfile in allfiles:

    if os.path.splitext(eachfile)[1] =='.tif':

        allhdffiles.append(eachfile)

#首先你需要新建一个Excel并建好几个工作表，当然你写一个改变之后的单元格参量让数据都写一个工作表里也行，但是我不推荐，容易错

 wb = openpyxl.load_workbook('D:/市NDVI.xlsx')#excel路径

 ws = wb["陕西"] #工作表名字

 ws.cell(row = 1, column = w+1, value = workfiles_province[w]) #我这里把一个省的几个市的数据写在同一个工作表里面，所以先写个头，row行数，column列数，value你要写入的值



 #主体循环，计算，一个文件计算一个值，写入一个单元格

 while n<len(allhdffiles):

  i=0

  k=0

  q=0

#打开文件

  dataset=gdal.Open(inpath+"/"+allhdffiles[n])

#栅格矩阵的列数

  im_width = dataset.RasterXSize

 #print(im_width)

#栅格矩阵的行数

  im_height = dataset.RasterYSize

 #print(im_height)

#读取某一像素点的值

#（1）读取一个波段，其参数为波段的索引号，波段索引号从1开始(我打开的这幅图像只有一个波段)

  band=dataset.GetRasterBand(1)

#（2）用ReadAsArray(<xoff>, <yoff>, <xsize>, <ysize>)，读出从(xoff,yoff)开始，大小为(xsize,ysize)的矩阵。以下为读取整幅图像

  im_datas=band.ReadAsArray(0,0,im_width,im_height)

  #两个循环，遍历所有像元，剔除Nodata，k为有效像元值和，q为有效像元个数

  while i<im_height:

   j=0

   while j<im_width:

    if(im_datas[i,j]!=-32768 and im_datas[i,j]!=65535):

      k+=im_datas[i,j]

#此处可能会给你报错，提示K值溢出什么的，把k=0改成k=0.00就可以
#（把k定义为浮点数，要是仍提示溢出，那就是你没有成功的把K定义成浮点数）
#要是你处理的数据都超过了浮点数的上限，我觉得你不会在百度上寻找解决办法（狗头）

      q+=1

    j+=1

   i+=1

  del dataset #加和完成后删除读入的文件，释放储存

  ws.cell(row = n+2, column = w+1, value = k/q) #写入单元格

  n+=1

 print(workfiles_province[w]) #当一个市完成后，输出一下市名，不然那么长时间啥也不输出我怎么知道这程序活着没有

 w+=1

 wb.save('D:/市NDVI.xlsx')#写完一个工作表保存一次。

在此感谢我数据处理期间在百度上浏览过的所有有用文章提供的零碎知识，参考较多的如下：

百度文库提供的《下载安装运行HEG处理modis数据,投影转换》

下载安装运行HEG处理modis数据,投影转换 - 百度文库

科学网李旭的《MODIS产品命名规则》

科学网—MODIS产品命名规则 - 李旭的博文

ASDN letskillthisgis的《ArcMap|批量按掩膜提取（模型构建器）》

ArcMap|批量按掩膜提取（模型构建器）_jilliranika-程序员宅基地

博客园岁时的《python调用HEG工具批量处理MODIS数据》

python调用HEG工具批量处理MODIS数据 - 岁时 - 博客园

菜鸟教程的编程基础

菜鸟教程 - 学的不仅是技术，更是梦想！