最近涉及到MODIS 植被指數(shù)、地表溫度、蒸散等時序數(shù)據(jù)產(chǎn)品的處理桶良，因此對MODIS數(shù)據(jù)批處理的方法進行了改進。主要是想充分利用工作站的性能沮翔，提升數(shù)據(jù)處理的效率陨帆。雖然用了兩三天時間，但是”磨刀不誤砍柴工“采蚀，還是值得的疲牵。
批處理工具MRT
MODIS數(shù)據(jù)批處理大多是依靠MRT工具，通過執(zhí)行cmd命令的方式完成榆鼠，需要提前手動生成相關(guān)產(chǎn)品處理的prm文件纲爸，可參考：
https://blog.csdn.net/jingchenlin1996/article/details/88548728
如果不同年份和月份的數(shù)據(jù)存放在不同的目錄下，那么就涉及到了文件夾的遍歷妆够。我們選擇了 python 的os.walk()函數(shù)识啦，并通過os.popen()函數(shù)執(zhí)行批處理的cmd命令。
#! usr/bin/python
#coding=utf-8
# 運行此文件前需要先運行MRT生成prm配置文件
import os
def MRTBacth(fileInpath,fileOutpath,prmFilepath):
for dirPath, dirNames, fileNames in os.walk(fileInpath):
if len(fileNames)> 0:
strCMD = 'cd .. && C: && cd C:/MRT/bin &&java -jar MRTBatch.jar -d ' + dirPath + \
' -p prmFilepath -o ' + fileOutpath + '&&MRTBatch.bat'
res = os.popen(strCMD)
print(res.read())
print(fileOutpath +': Finished')
MRTBacth('H:/DATA/MOD16A2/2001')
到這里都只是簡單的單核計算神妹。
并行計算
python并行計算可參考資料：
以及基礎(chǔ)的并發(fā)颓哮、并行、阻塞等概念鸵荠，參考：
Python+MRT 并行處理
遇見的問題
我們選擇了multiprocessing的多進程包冕茅，在上面MRTBatch()函數(shù)的基礎(chǔ)上增加下面的內(nèi)容：
#導(dǎo)入multiprocessing包和functools包的partial模塊
import multiprocessing as mp
from functools import partial
#windows環(huán)境下主程序應(yīng)放在if __name__ =="__main__"下
if __name__ == "__main__":
fileInpath = 'H:/DATA/MOD16A2/'
fileOutpath='H:\DATA\RESULT\MOD16A2'
prmFile='H:/DATA/MOD16A2.prm'
#獲取輸入路徑下的所有包含文件的文件路徑
fileList=GetFile(fileInpath)
#創(chuàng)建10個進程的進程池，
pool = mp.Pool(10)
#因為MRTBactch()需要傳入多參數(shù)蛹找，可通過partial()函數(shù)固定其他參數(shù)
func=partial(MRTBacth,fileOutpath=fileOutpath,prmFile=prmFile)
#通過pool.map()函數(shù)實現(xiàn)并行計算
res=pool.map(func,fileList)
print(res.get())
pool.close()
pool.join()
看似沒什么問題姨伤，實際上計算的結(jié)果并不對。正常輸出的tif文件數(shù)據(jù)量約180M熄赡，但此時程序算出來的數(shù)據(jù)量約1.7G姜挺，差不多10輩的差距。更意外的是只生成了一個TmpMosaic.hdf 的中間過程文件彼硫。
通過查看日志文件凌箕，發(fā)現(xiàn)雖然進程池中的進程都執(zhí)行了mrtmosaic.exe和resample.exe，但是這倆exe在不同進程中的輸入輸出路徑都是一樣的词渤。
后來發(fā)現(xiàn)當(dāng)執(zhí)行MRT批處理中“java -jar MRTBatch.jar -d ' + dirPath + -p prmFilepath -o ' + fileOutpath ”的命令時，MRT會在其安裝目錄下生成文件名為MRTBatch.bat的文件缺虐，文件內(nèi)容如下（請忽略具體的輸出輸出路徑）：
mrtmosaic -s "1 0 0 0 0" -i H:\DATA\RESULT\MOD16A2/2001\353\MOD16A2.A2001353_mosaic.prm -o H:\DATA\RESULT\MOD16A2/2001\353\TmpMosaic.hdf
resample -p H:\DATA\RESULT\MOD16A2/2001\353\MOD16A2.A2001353_resample.prm
del H:\DATA\RESULT\MOD16A2/2001\353\TmpMosaic.hdf
原來MRTBatch.bat的文件記錄的時后續(xù)調(diào)用mrtmosaic.exe和resample.exe所需的輸入?yún)?shù)芜壁。這解釋了為什么MRT批處理命令中要有”MRTBatch.bat“這一命令，這也解釋了為什么不同進程的輸出結(jié)果是相同的高氮。因為每個進程都會在同一目錄下生成MRTBatch.bat文件慧妄，文件會自動覆蓋，最終只記錄了最后一次的內(nèi)容剪芍。
解決辦法
改進的方法就是獨立調(diào)用mrtmosaic.exe和resample.exe塞淹。通過查看MRT幫助文檔，了解了兩個exe的參數(shù)罪裹，對程序進行了調(diào)試和修改饱普。最終的程序代碼如下：
#! usr/bin/python
#coding=utf-8
# 運行此文件前需要先運行MRT生成prm配置文件
import os,shutil,time,re
import multiprocessing as mp
from functools import partial
#文件夾遍歷，獲取包含hdf文件的所有文件夾状共，并以字符串列表返回
def GetFile(fileInpath):
fileList=[]
for dirPath, dirNames, fileNames in os.walk(fileInpath):
if len(fileNames) > 0:
fileList.append(dirPath)
return fileList
#獲取某文件夾下所有后綴為hdf的文件名套耕，以字符串列表返回
def GetFileName(fileInpath,fileType):
fileNames = os.listdir(fileInpath)
fileNameArr = []
for fileName in fileNames:
if os.path.splitext(fileName)[1] == fileType:
fileNameArr.append(fileName)
return fileNameArr
#mrtmosaic.exe
def MrtmosaicPro(fileInpath,fileOutpath,SpectralSubset):
#build mosaic prm file，which records all the datafiles under the fileInpath
fileNames = GetFileName(fileInpath, '.hdf')
prmName = fileOutpath + '/' + 'tmpMosaic.prm'
for fileName in fileNames:
with open(prmName, 'a+') as f:
f.write(fileInpath+'/'+fileName + '\n')
#get the product type and DOY information
productType=fileNames[0].split('.')[0]
DOY=fileNames[0].split('.')[1]
fileOutName=fileOutpath +'/'+productType+'.'+DOY+'.'+'TmpMosaic.hdf'
#run mrtmosaic.exe
strCMD = 'cd .. && C: && cd C:/MRT/bin && mrtmosaic -i '+prmName+ ' -s "'+SpectralSubset+'" -o '+fileOutName
res = os.popen(strCMD)
print(res.read())
print(fileOutName +': Finished')
return fileOutName,prmName
#fileInpath:MOD16A2.A2010001.tmpmosaic.hdf
#fileOutpath:
#prmFilepath:
#SpectralSubset:'1 0 0 0 0'
def ResamplePro(fileInpath,fileOutpath,prmFilepath,SpectralSubset):
#copy the prm file into fileoutpath and update it
shutil.copy(prmFilepath,fileOutpath)
filedata=[]
newPrmFilepath=fileOutpath+'/'+os.path.basename(prmFilepath)
fileOutName=os.path.basename(fileInpath).split('.')[0]+'.'+os.path.basename(fileInpath).split('.')[1]+'.tif'
targetStr1='INPUT_FILENAME = '
newStr1='INPUT_FILENAME = '+fileInpath
targetStr2='SPECTRAL_SUBSET = '
newStr2 = 'SPECTRAL_SUBSET = (' + str(SpectralSubset.count('1'))+')'
targetStr3 = 'OUTPUT_FILENAME = '
newStr3 = 'OUTPUT_FILENAME = ' + fileOutpath+'/'+fileOutName
#read in the prmfile and replace the variables of input filename,SPECTRAL_SUBSET and OUTPUT_FILENAME
with open(newPrmFilepath, "r") as file:
#lines=file.read()
for line in file:
if re.match(targetStr1,line):
line=newStr1
if re.match(targetStr2,line):
line=newStr2
if re.match(targetStr3,line):
line=newStr3
filedata.append(line)
with open(newPrmFilepath,'a+') as file:
file.seek(0)
file.truncate()
for i in range(len(filedata)):
file.write(filedata[i]+'\n')
#run resample.exe
strCMD = 'cd .. && C: && cd C:/MRT/bin && resample -p ' +newPrmFilepath
res = os.popen(strCMD)
print(res.read())
print(fileOutName + ': Finished')
return newPrmFilepath
def MRTBacth(fileInpath,fileOutpath,prmFile,SpectralSubset):
#根據(jù)實際情況自己改
newFileOutpath=fileOutpath+fileInpath[15:24]
if os.path.exists(newFileOutpath) == False:
os.makedirs(newFileOutpath)
result= MrtmosaicPro(fileInpath, newFileOutpath, SpectralSubset)
#use the output from mrtmosaic.exe as the inputfile to resample.exe
newPrmFilepath=ResamplePro(result[0], newFileOutpath, prmFile, SpectralSubset)
#delete prm file and tmpmosaic.hdf
strDelCMD='del ' + result[0] + '&& del ' + result[1]+'&& del '+newPrmFilepath
res = os.popen(strDelCMD)
if __name__ == "__main__":
startTime = time.time()
fileInpath = 'H:/DATA/MOD16A2/'
fileOutpath='H:\DATA\RESULT\MOD16A2'
prmFile='H:/DATA/MOD16A2.prm'
SpectralSubset="1 0 0 0 0"
fileList=GetFile(fileInpath)
#MRTBacth(fileList[0],fileOutpath,prmFile,SpectralSubset)
# cores number
cores = mp.cpu_count()
pool = mp.Pool(int(cores * 2 / 3))
func=partial(MRTBacth,fileOutpath=fileOutpath,prmFile=prmFile,SpectralSubset=SpectralSubset)
res=pool.map(func,fileList)
print(res.get())
pool.close()
pool.join()
endTime = time.time()
print("used time is ", endTime - startTime)
總結(jié)
以前使用MRT沒有注意到程序運行的一些細(xì)節(jié)峡继，通過這次實驗冯袍，了解了MRT內(nèi)部運行的機制。機器的內(nèi)存和顯卡資源還未利用起來鬓椭，后續(xù)將對程序進一步優(yōu)化和改進颠猴。最后放一張cpu的圖：
