Python GDAL批处理MODIS产品

最新推荐文章于 2023-04-03 10:43:43 发布

江北20190411

最新推荐文章于 2023-04-03 10:43:43 发布

阅读量2.6k

点赞数 10

分类专栏： python 遥感

本文链接：https://blog.csdn.net/tk20190411/article/details/105975593

版权

遥感同时被 2 个专栏收录

40 篇文章 33 订阅

订阅专栏

python

36 篇文章 42 订阅

订阅专栏

前言

学习Python GDAL有几天了,这次把常用的功能写成批量处理的函数，方便以后可以直接使用。
事实上,前面的几何校正更推荐使用MCTK和MRT进行处理,毕竟有去除bowtie的算法。
在这里使用的数据是MOD13A3全球数据,该数据提供每月1km分辨率的3级正弦曲线投影网格产品。
之后打算实现将12副图绘制在一个图窗中。

完整代码

from osgeo import gdal, osr
from osgeo import gdal_array
import os
import time
import numpy as np

start = time.time() 
'''
脚本目的:批量读取tif文件并进行几何校正,定标和投影转换,裁剪
'''
'''
第一个函数:批量读取文件夹中hdf格式的modis产品数据并进行几何校正(WGS84)
'''
def readHdfWithGeo(hdfFloder, saveFloder):
    #获取输入文件夹中所有的文件名
    hdfNameList = os.listdir(hdfFloder)
    #遍历文件名列表中所有文件
    for i in range(len(hdfNameList)):
        #获取文件名后缀
        dirname, basename = os.path.split(hdfNameList[i])
        filename, txt = os.path.splitext(basename)
        #判断文件后缀是否为 .hdf
        if txt == '.hdf':
            hdfPath = hdfFloder + os.sep + hdfNameList[i]
            #打开hdf文件
            datasets = gdal.Open(hdfPath)
            #打开子数据集
            dsSubDatasets = datasets.GetSubDatasets()
            #打开ndvi数据
            ndviRaster = gdal.Open(dsSubDatasets[0][0])
            #获取元数据
            metaData = datasets.GetMetadata()
            #for key, value in metaData.items():
            #   print('{key}:{value}'.format(key = key, value = value))
            #获取数据时间
            time = metaData['RANGEBEGINNINGDATE']
            #命名输出完整路径文件名
            outName = saveFloder + os.sep + time + '.tif'
            #进行几何校正
            geoData = gdal.Warp(outName, ndviRaster,
                                dstSRS = 'EPSG:4326', format = 'GTiff',
                                resampleAlg = gdal.GRA_Bilinear)
            del geoData
            print('{outname} deal end'.format(outname = outName))
            
start = time.time()      
    
readHdfWithGeo(r'F:\GDAL learning\MODIS\MOD13A3\2012', 
               r'F:\GDAL learning\MODIS\MOD13A3\NDVItif_withGeo')
end = time.time()
print('deal spend: {s} s'.format(s = end-start))


'''
第二个函数:批量读取tif文件(适用于多波段tif文件)存入数组,并进行定标
'''
def readTifAsArray(tifPath):
    dataset = gdal.Open(tifPath)     
    if dataset == None:
        print(tifPath + "文件错误")
        return tifPath
            
    image_datatype = dataset.GetRasterBand(1).DataType  
    row = dataset.RasterYSize
    col = dataset.RasterXSize
    nb  = dataset.RasterCount
    proj = dataset.GetProjection()
    gt = dataset.GetGeoTransform()
    
    if nb != 1:
        array = np.zeros((row, col, nb), 
                        dtype = gdal_array.GDALTypeCodeToNumericTypeCode(
                                image_datatype))
        for b in range(nb):
            band = dataset.GetRasterBand(b + 1)
            nan = band.GetNoDataValue()
            array[:, :, b] = band.ReadAsArray()
    else:
        array = np.zeros((row,col),
                         dtype = gdal_array.GDALTypeCodeToNumericTypeCode(
                         image_datatype))
        band = dataset.GetRasterBand(1)        
        nan = band.GetNoDataValue()        
        array = band.ReadAsArray()
    return array, nan, gt, proj

'''  
第三个函数:写出tif文件
'''
def writeTiff(im_data, nan, im_geotrans, im_proj, path):
    if 'int8' in im_data.dtype.name:
        datatype = gdal.GDT_Byte
    elif 'int16' in im_data.dtype.name:
        datatype = gdal.GDT_UInt16
    else:
        datatype = gdal.GDT_Float32
    if len(im_data.shape) == 3:
        im_bands, im_height, im_width = im_data.shape
    elif len(im_data.shape) == 2:
        im_data = np.array([im_data])
        im_bands, im_height, im_width = im_data.shape
   
    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(path, int(im_width), int(im_height), int(im_bands), datatype)
    
    if(dataset != None):
        dataset.SetGeoTransform(im_geotrans) #写入仿射变换参数
        dataset.SetProjection(im_proj) #写入投影
    for i in range(im_bands):
        dataset.GetRasterBand(i+1).WriteArray(im_data[i])
        outBand = dataset.GetRasterBand(i + 1)
        outBand.FlushCache()
        outBand.SetNoDataValue(nan)
    del dataset
 
'''
第四个函数:批量进行定标计算
'''
def cal(tifFloder, saveFloder, scale):
    tifNameList = os.listdir(tifFloder)  

    for i in range(len(tifNameList)):
        
        filename, txt = os.path.splitext(tifNameList[i])
        
        if txt == '.tif':
            tifPath = tifFloder + os.sep +  tifNameList[i]
            array = readTifAsArray(tifPath)
            outName = saveFloder + os.sep + filename + '_cal' + txt
            calData = writeTiff(array[0]/scale, array[1]/scale, 
                            array[2], array[3], outName)
            print('{filename} deal end '.format(filename = outName))
            del calData
            
start = time.time() 
cal(r'F:\GDAL learning\MODIS\MOD13A3\NDVItif_withGeo',
    r'F:\GDAL learning\MODIS\MOD13A3\NDVItif_cal', 10000)        
end = time.time()
print('deal spend: {s} s'.format(s = end-start))

'''
第五个函数:批量进行重投影
'''    
def reproject(tifFloder, saveFloder, proj4):
    tifNameList = os.listdir(tifFloder)  

    srs = osr.SpatialReference()
    srs.ImportFromProj4(proj4)

    for i in range(len(tifNameList)):

        filename, txt = os.path.splitext(tifNameList[i])
        
        if txt == '.tif': 
            tifPath = tifFloder + os.sep +  tifNameList[i]
            outName = saveFloder + os.sep + filename + '_rep' + txt
            reproData = gdal.Warp(outName, tifPath, dstSRS = srs, 
                                  xRes = 1000, yRes = 1000,          
                                  resampleAlg = gdal.GRA_Bilinear,
                                  outputType=gdal.GDT_Float32)
            print('{filename} deal end '.format(filename = outName))
            del reproData
            
start = time.time()       
reproject(r'F:\GDAL learning\MODIS\MOD13A3\NDVItif_cal',
          r'F:\GDAL learning\MODIS\MOD13A3\NDVItif_albers',
          '+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs')
end = time.time()
print('deal spend: {s} s'.format(s = end-start))

'''
第六个函数:批量进行裁剪
'''
def cut(tifFloder, saveFloder, shpFile):
    tifNameList = os.listdir(tifFloder)  

    for i in range(len(tifNameList)):

        filename, txt = os.path.splitext(tifNameList[i])
        
        if txt == '.tif': 
            tifPath = tifFloder + os.sep +  tifNameList[i]
            outName = saveFloder + os.sep + filename + '_cut' + txt
            cut_ds = gdal.Warp(outName, tifPath,
                   cutlineDSName = shpFile,
                   cropToCutline = True)
            print('{filename} deal end '.format(filename = outName))
            del cut_ds
 
start = time.time()          
cut(r'F:\GDAL learning\MODIS\MOD13A3\NDVItif_albers',
        r'F:\GDAL learning\MODIS\MOD13A3\NDVItif_wrap',
        r'F:\GDAL learning\shp\BeiJingBei.shp')
end = time.time()
print('deal spend: {s} s'.format(s = end-start))

总共用了20+s

F:\GDAL learning\MODIS\MOD13A3\NDVItif_withGeo\2012-01-01.tif deal end
F:\GDAL learning\MODIS\MOD13A3\NDVItif_withGeo\2012-02-01.tif deal end
F:\GDAL learning\MODIS\MOD13A3\NDVItif_withGeo\2012-03-01.tif deal end
F:\GDAL learning\MODIS\MOD13A3\NDVItif_withGeo\2012-04-01.tif deal end
F:\GDAL learning\MODIS\MOD13A3\NDVItif_withGeo\2012-05-01.tif deal end
F:\GDAL learning\MODIS\MOD13A3\NDVItif_withGeo\2012-06-01.tif deal end
F:\GDAL learning\MODIS\MOD13A3\NDVItif_withGeo\2012-07-01.tif deal end
F:\GDAL learning\MODIS\MOD13A3\NDVItif_withGeo\2012-08-01.tif deal end
F:\GDAL learning\MODIS\MOD13A3\NDVItif_withGeo\2012-09-01.tif deal end
F:\GDAL learning\MODIS\MOD13A3\NDVItif_withGeo\2012-10-01.tif deal end
F:\GDAL learning\MODIS\MOD13A3\NDVItif_withGeo\2012-11-01.tif deal end
F:\GDAL learning\MODIS\MOD13A3\NDVItif_withGeo\2012-12-01.tif deal end
deal spend: 6.361145496368408 s
F:\GDAL learning\MODIS\MOD13A3\NDVItif_cal\2012-01-01_cal.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_cal\2012-02-01_cal.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_cal\2012-03-01_cal.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_cal\2012-04-01_cal.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_cal\2012-05-01_cal.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_cal\2012-06-01_cal.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_cal\2012-07-01_cal.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_cal\2012-08-01_cal.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_cal\2012-09-01_cal.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_cal\2012-10-01_cal.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_cal\2012-11-01_cal.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_cal\2012-12-01_cal.tif deal end 
deal spend: 3.9119620323181152 s
F:\GDAL learning\MODIS\MOD13A3\NDVItif_albers\2012-01-01_cal_rep.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_albers\2012-02-01_cal_rep.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_albers\2012-03-01_cal_rep.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_albers\2012-04-01_cal_rep.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_albers\2012-05-01_cal_rep.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_albers\2012-06-01_cal_rep.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_albers\2012-07-01_cal_rep.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_albers\2012-08-01_cal_rep.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_albers\2012-09-01_cal_rep.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_albers\2012-10-01_cal_rep.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_albers\2012-11-01_cal_rep.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_albers\2012-12-01_cal_rep.tif deal end 
deal spend: 8.790398836135864 s
F:\GDAL learning\MODIS\MOD13A3\NDVItif_wrap\2012-01-01_cal_rep_cut.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_wrap\2012-02-01_cal_rep_cut.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_wrap\2012-03-01_cal_rep_cut.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_wrap\2012-04-01_cal_rep_cut.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_wrap\2012-05-01_cal_rep_cut.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_wrap\2012-06-01_cal_rep_cut.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_wrap\2012-07-01_cal_rep_cut.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_wrap\2012-08-01_cal_rep_cut.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_wrap\2012-09-01_cal_rep_cut.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_wrap\2012-10-01_cal_rep_cut.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_wrap\2012-11-01_cal_rep_cut.tif deal end 
F:\GDAL learning\MODIS\MOD13A3\NDVItif_wrap\2012-12-01_cal_rep_cut.tif deal end 
deal spend: 1.5061616897583008 s

江北20190411

关注

10
点赞
踩
52

收藏

觉得还不错? 一键收藏
2
评论
Python GDAL批处理MODIS产品

批量几何校正MODIS产品from osgeo import gdal, osrimport osimport time'''批量读取文件夹中hdf格式的modis产品数据并进行几何校正(WGS84)'''def readHdfWithGeo(hdfFloder, saveFloder): #获取输入文件夹中所有的文件名 hdfNameList = os.listdi...
复制链接

扫一扫