python調用HEG工具批量處理MODIS數據

下麵的代碼主要用於使用python語言調用NASA官方的MODIS處理工具HEG進行投影坐標轉換與重採樣批量處理
主要參考

1. HEG的用戶手冊：https://newsroom.gsfc.nasa.gov/sdptoolkit/HEG/HEG215/EED2-TP-030_Rev01_HEG_UsersGuide_2.15.pdf
2. HEG批處理幫助：https://newsroom.gsfc.nasa.gov/sdptoolkit/HEG/HEG_Batch_job_Help.htm

主要的註意事項如下：

1. 首先按照官方指南安裝HEG工具，具體安裝步驟可參考我的上篇博客：https://www.cnblogs.com/yhpan/p/11298595.html
2. 根據HEG用戶手冊批量生成批處理參數文件，可以在HEG工具中生成一個文件，拿來自己改改用用
3. 具體調用哪一個工具，參數文件如何寫，請一定仔細閱讀用戶手冊，東西全都在上面。一般常用的是resample.exe和swtif.exe，如果實在無法判斷可以先用HEG的GUI處理一個自己的數據，保存一個prm文件，然後根據這個文件中的參數，對照著用戶手冊一個一個的找，就可以了。
4. 生成參數文件寫入時一定要註意，設定換行符為‘\n’,fo=open(prmfilename,'w',newline='\n')，否則由於在windows系統下預設換行符為‘\r\n’,程式無法運行成功

下麵是源碼分享

# -*- coding: utf-8 -*-
"""
Created on Sun Feb 16 11:27:19 2020
調用HEG相關工具批處理MODIS數據，主要完成投影坐標轉換與重採樣
@author: pan
"""
import os

# 設置HEG相關環境變數
os.environ['MRTDATADIR']='D:/MyApps/HEG/HEG_Win/data'
os.environ['PGSHOME']='D:/MyApps/HEG/HEG_Win/TOOLKIT_MTD'
os.environ['MRTBINDIR']='D:/MyApps/HEG/HEG_Win/bin'

# 設置HEG的bin路徑
hegpath = 'D:/MyApps/HEG/HEG_Win/bin'
# 指定處理模塊的可執行程式文件路徑，此處採用resample.exe，可以根據具體的處理問題設置
hegdo = os.path.join(hegpath, 'resample.exe')
hegdo = hegdo.replace('\\', '/') # 全路徑以“/”連接

# 指定輸入數據的路徑
inpath = r'C:\Users\pan\Desktop\Py_ex\data\hdf'
inpath = inpath.replace('\\', '/')
# 指定輸出數據的路徑
outpath = r'C:\Users\pan\Desktop\Py_ex\data\hdf\out'
outpath = outpath.replace('\\', '/')
# os.chdir(inpath) #改變當前工作目錄到輸入數據目錄

# 獲取當前文件夾下的所有hdf文件
allfiles = os.listdir(inpath)
allhdffiles = []
for eachfile in allfiles:
    if os.path.splitext(eachfile)[1] =='.hdf':
        allhdffiles.append(eachfile)
print('--'*20)
print('文件數量為:', len(allhdffiles),',所有hdf文件如下')
print('  '+'\n  '.join(allhdffiles))
print('--'*20)

# prm文件設置模塊，需要首先在HEG工具中生成一個參考的prm文件，示例如下
# 設置prm文件存儲路徑
prmpath = r"C:\Users\pan\Desktop\Py_ex\data\hdf\prm"
prmpath = prmpath.replace('\\', '/')
for eachhdf in allhdffiles:
    prm=['NUM_RUNS = 1\n',
      'BEGIN\n',
      'INPUT_FILENAME = ' + inpath+'/'+eachhdf+'\n',
      'OBJECT_NAME = MODIS_Grid_8Day_1km_LST|\n',
      'FIELD_NAME = LST_Day_1km\n',
      'BAND_NUMBER = 1\n','SPATIAL_SUBSET_UL_CORNER = ( 90.0 -180.0 )\n',
      'SPATIAL_SUBSET_LR_CORNER = ( -90.0 180 )\n',
      'RESAMPLING_TYPE = BI\n',
      'OUTPUT_PROJECTION_TYPE = ALBERS\n',
      'ELLIPSOID_CODE = WGS84\n',
      'OUTPUT_PROJECTION_PARAMETERS = ( 0.0 0.0 25.0 47.0 105.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0  )\n',
      'OUTPUT_PIXEL_SIZE = 500.0\n',
      'OUTPUT_FILENAME = ' + outpath+'/'+eachhdf+'_out.tif\n',
      'OUTPUT_TYPE = GEO\n',
      'END\n']
    prmfilename=prmpath +'/'+ eachhdf+'.prm'
    prmfilename=prmfilename.replace('\\', '/')
    #這裡一定要註意，設定換行符為‘\n’,否則由於在windows系統下預設換行符為‘\r\n’,則無法運行成功
    fo=open(prmfilename,'w',newline='\n')
    fo.writelines(prm)
    fo.close()

for eachhdf in allhdffiles:
    prmfilepath=prmpath +'\\'+ eachhdf + '.prm'
    prmfilepath=prmfilepath.replace('\\', '/')
    try:
        resamplefiles = '{0} -P {1}'.format(hegdo, prmfilepath)
        os.system(resamplefiles)        
        print(eachhdf + ' has finished')
    except:
        # 提示錯誤信息
        print(eachhdf + 'was wrong')