Maiac文件的AOD信息提取-by python

Modis数据向来风骚，其HDF里包含了很多subdataset，其中有一个maiac的数据尤为特别。Maiac文件里含有大概12个数据集，每个数据集里又有4个波段（维度），如果按GDAL的translate函数直接转换，将得到错误的结果，会只得到第一个波段Band1如下：

--------------------------------------------------------分割线------------

导致其他波段无法表达，这就不合适了。所以就需要重新理清一下顺序。

首先安装GDAL，具体教程可以百度，但是有个注意的是安装时请使用typical模式，不要complete，否则会出错。

接着使用GDAL的translate函数转换出一张影像图：

import subprocess

subprocess.call('gdal_translate'+' -sds'+' -b'+' 2'+' D:/Thesis/ML/aod/MOD13Q1.A2017081.h28v06.006.2017111085049.hdf'+" D:\Thesis\ML\aodband2\maiac.tif")

其中（sds）这些变量前面要加空格，否则语法不对~

好了，这下可以开始处理maiac了，思路是把4个波段的最大值都赋值到一张图上，最大限度的利用数据。。。。

import gdal, osr

import numpy as np

ds = gdal.Open('D:/Thesis/ML/aod/MCD19A2.A2018001.h28v06.006.2018121012322.hdf')

sub=ds.GetSubDatasets()

aod55=gdal.Open(sub[1][0]).ReadAsArray()

a0=aod55[0,:]

a1=aod55[1,:]

a2=aod55[2,:]

a3=aod55[3,:]

[cols, rows] = a0.shape

a4=np.zeros((cols,rows))

for i in range(len(a0)):

for j in range(len(a0[1])):

a4[i,j]=max(a0[i,j],a1[i,j],a2[i,j],a3[i,j])

ds2 = gdal.Open("D:\Thesis\ML\maiac_02.tif")

band = ds2.GetRasterBand(1)

arr = band.ReadAsArray()

[cols, rows] = arr.shape

format = "GTiff"

driver = gdal.GetDriverByName(format)

outDataRaster = driver.Create("D:/Thesis/ML/aodband2/aod4.tif", rows, cols, 1, gdal.GDT_Int16)

outDataRaster.SetGeoTransform(ds2.GetGeoTransform())##sets same geotransform as input

outDataRaster.SetProjection(ds2.GetProjection())##sets same projection as input

outDataRaster.GetRasterBand(1).WriteArray(a4)

outDataRaster.FlushCache()

del outDataRaster

最后得到的是：

跟第一张图还是有点差别的。。。因为有数据填补。