MOD043km的hdf数据转为tif

MOD043km的hdf数据转为tifc尝试了两种办法，1是靠

gdal.open(r'D:/Thesis/ML/modis3km/MOD04_3K.A2018001.0320.061.2018003202214.hdf')

代码如下：

import gdal, osr

import numpy as np

import os

import rasterio

from rasterio.warp import calculate_default_transform, reproject, Resampling

from rasterio import crs

from pymodis.convertmodis_gdal import convertModisGDAL

ds = gdal.Open(r'D:/Thesis/ML/modis3km/MOD04_3K.A2018001.0320.061.2018003202214.hdf')

sub=ds.GetSubDatasets()

aod=gdal.Open(sub[10][0]).ReadAsArray()##这步有个问题aod=gdal.Open(sub[10]).ReadAsArray()的时候立刻报错，还是要[10][0]才可以，后续的问题不知道是不是这步有毛病

band = ds.GetRasterBand(1)#2

#arr = band.ReadAsArray()#3，

#[cols, rows] = arr.shape

[cols, rows] = aod.shape

format = "GTiff"#5

driver = gdal.GetDriverByName(format)#6

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

outDataRaster.SetGeoTransform(ds.GetGeoTransform())##sets same geotransform as input8

outDataRaster.SetProjection(ds.GetProjection())##sets same projection as input9

outDataRaster.GetRasterBand(1).WriteArray(aod)#10，5-10都是赋值给新数据

outDataRaster.FlushCache()#把缓存清理，很重要

del outDataRaster

结果：

是有问题的，它按行列号重新显示了图像，没有坐标系，

ds.GetGeoTransform()也有问题，其中6个参数分别为：x0, dx, dxdy, y0, dydx, dy = ds.GetGeoTransform()，有分辨率，坐标参数，还有行列数目貌似

于是尝试第二种办法：

import subprocess

import os

subprocess.call('gdal_translate'+' -of GTiff'+' "HDF4_EOS:EOS_SWATH:"D:/Thesis/ML/modis3km/MOD04_3K.A2018001.0320.061.2018003202214.hdf":mod04:Image_Optical_Depth_Land_And_Ocean"'+''+' "D:/Thesis/ML/aodband2/aod5.tif"')

一定要注意其中的空格，还有文件名前面的引号

结果：

这个结果可以接受了，是想要的wgs84