利用gdal、rasterio将modis文件进行格式转换、投影转换

modis的hdf文件在存储上有优势，但是在实际使用过程中存在一定的弊端。例如本次重点讲解的NDVI-16D-1km、MAIAC-1D-1km两类文件，其中maiac的文件在aod属性中包含4个波段的文件，需要将4个波段最大化的利用起来。而ndvi文件则较为简单只需要把hdf文件中对应部分取用即可。

1、MAIAC文件，在前文已经有所讲解：需要注意的是，要先备一份带有坐标系的tif文件，其实如果对gdal很熟悉的话，这部分是可以跳过的。所以本次的代码任然有优化和改进的空间，但是感觉在hdf转tif这部中rasterio的效率比gdal高多了

import gdal, osr

import numpy as np

import os

import rasterio

from rasterio.warp import calculate_default_transform, reproject, Resampling

from rasterio import crs

path= r'D:/Thesis/fujian/201711-201806/maiac/rawdata'

ou=r'D:/Thesis/fujian/201711-201806/maiac/tif'

filers=os.listdir(path)

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

band = ds2.GetRasterBand(1)#2

arr = band.ReadAsArray()#3，

[cols, rows] = arr.shape#4，1234里都是把预设好那个gdal导出的数据做参数提取

format = "GTiff"#5

driver = gdal.GetDriverByName(format)#6

for u in filers:

inpute=path+'/'+u

outp=ou+'/'+u+'.tif'

ds=gdal.Open(inpute)

sub=ds.GetSubDatasets()

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

if aod55.ndim == 4:

a0=aod55[0,:]#提取第一个波段（维度）

a1=aod55[1,:]#2

a2=aod55[2,:]#3

a3=aod55[3,:]#4

[col, row] = a0.shape

a4=np.zeros((col,row))#建立一个空的数组，到时候放新的数据

for i in range(len(a0)):###只能用range和长度来循环i和j，本来想说通过索引。。。不知道为什么不行，导致你以下的数据必须一致才能继续用

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

a4[i,j]=max(a0[i,j],a1[i,j],a2[i,j],a3[i,j],a4[i,j])##挑里面最大的，如果还有a4的话，那最小值就是0了

outDataRaster = driver.Create(outp, rows, cols, 1, gdal.GDT_Int16)#7gdal.GDT_Int16还有其他参数，比如gdal.GDT_Byte是8bit

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

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

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

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

del outDataRaster

elif aod55.ndim == 3:

a0=aod55[0,:]#提取第一个波段（维度）

a1=aod55[1,:]#2

a2=aod55[2,:]#3

# a3=aod55[3,:]#4

[col, row] = a0.shape

a4=np.zeros((col,row))#建立一个空的数组，到时候放新的数据

for i in range(len(a0)):###只能用range和长度来循环i和j，本来想说通过索引。。。不知道为什么不行，导致你以下的数据必须一致才能继续用

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

a4[i,j]=max(a0[i,j],a1[i,j],a2[i,j],a4[i,j])##挑里面最大的，如果还有a4的话，那最小值就是0了

outDataRaster = driver.Create(outp, rows, cols, 1, gdal.GDT_Int16)#7gdal.GDT_Int16还有其他参数，比如gdal.GDT_Byte是8bit

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

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

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

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

del outDataRaster

#投影转换，这部分既可以用rasterio也可以用subprocess.call调用gdal

repro = r'D:/Thesis/fujian/201711-201806/maiac/repro'

oufile=os.listdir(ou)

dst_crs = crs.CRS.from_epsg('32650')#4326是wgs84

for k in oufile:

src_img=ou+'/'+k

dst_img=repro+'/'+k

with rasterio.open(src_img) as src_ds:

profile = src_ds.profile

dst_transform, dst_width, dst_height = calculate_default_transform(src_ds.crs, dst_crs, src_ds.width, src_ds.height, *src_ds.bounds)

profile.update({'crs': dst_crs,'transform': dst_transform,'width': dst_width,'height': dst_height,'nodata': 0})

with rasterio.open(dst_img, 'w', **profile) as dst_ds:

for i in range(1, src_ds.count + 1):

src_array = src_ds.read(i)

dst_array = np.empty((dst_height, dst_width), dtype=profile['dtype'])

reproject(source=src_array,src_crs=src_ds.crs,src_transform=src_ds.transform,destination=dst_array,dst_transform=dst_transform,dst_crs=dst_crs,resampling=Resampling.cubic,num_threads=2)

dst_ds.write(dst_array, i)

2NDVI：

import gdal, osr

import numpy as np

import os

import rasterio

from rasterio.warp import calculate_default_transform, reproject, Resampling

from rasterio import crs

path= r'D:/minxinan/data/NDVI/ndvihdf'

ou=r"D:/minxinan/data/NDVI/ndvitif"

filers=os.listdir(path)

ds2 = gdal.Open("D:/minxinan/data/temp/ndvitemp.tif")#1

band = ds2.GetRasterBand(1)#2

arr = band.ReadAsArray()#3，

[cols, rows] = arr.shape#4，1234里都是把预设好那个gdal导出的数据做参数提取

format = "GTiff"#5

driver = gdal.GetDriverByName(format)#6

for u in filers:

inpute=path+'/'+u

outp=ou+'/'+u+'.tif'

ds=gdal.Open(inpute)

sub=ds.GetSubDatasets()

ndvi=gdal.Open(sub[0][0]).ReadAsArray()

outDataRaster = driver.Create(outp, rows, cols, 1, gdal.GDT_Int16)#7gdal.GDT_Int16还有其他参数，比如gdal.GDT_Byte是8bit

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

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

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

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

del outDataRaster

repro = r'D:/minxinan/data/NDVI/ndvirepro'

oufile=os.listdir(ou)

dst_crs = crs.CRS.from_epsg('4326')

for k in oufile:

src_img=ou+'/'+k

dst_img=repro+'/'+k

with rasterio.open(src_img) as src_ds:

profile = src_ds.profile

dst_transform, dst_width, dst_height = calculate_default_transform(src_ds.crs, dst_crs, src_ds.width, src_ds.height, *src_ds.bounds)

profile.update({'crs': dst_crs,'transform': dst_transform,'width': dst_width,'height': dst_height,'nodata': 0})

with rasterio.open(dst_img, 'w', **profile) as dst_ds:

for i in range(1, src_ds.count + 1):

src_array = src_ds.read(i)

dst_array = np.empty((dst_height, dst_width), dtype=profile['dtype'])

reproject(source=src_array,src_crs=src_ds.crs,src_transform=src_ds.transform,destination=dst_array,dst_transform=dst_transform,dst_crs=dst_crs,resampling=Resampling.cubic,num_threads=2)

dst_ds.write(dst_array, i)

import gdal, osr

import numpy as np

import os

filers=os.listdir(r'D:/minxinan/data/NDVI/ndvihdf')

for i in filers:

ds = gdal.Open('D:/minxinan/data/NDVI/ndvihdf/'+i)

sub=ds.GetSubDatasets()

aod=gdal.Open(sub[0][0]).ReadAsArray()

band = ds.GetRasterBand(1)#2

[cols, rows] = aod.shape

format = "GTiff"#5

driver = gdal.GetDriverByName(format)#6

outDataRaster = driver.Create(r"D:/minxinan/data/NDVI/ndvitif/"+i+'.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