遥感影像处理

学习遥感领域的小伙伴，经常会遇到处理多景影像的叠加运算，日均值合成，但是由于日数据像元存在Nodata空值，在计算均值过程中，又不想受到Nodata计算的影响。在学习数据处理的过程中不断学习，尝试用GDAL自我学习和记录笔记

# 导入相关的库
import os
import pandas as pd
import numpy as np
from osgeo import gdal

 #读取文件夹路径下影像       
def read_tif_file(path):
    filelist = os.listdir(path) 
    total_num = len(filelist)  
    
    image_array1 = np.empty((1200, 1200))
    for item in filelist:
        if item.endswith('.tif'):
            dataset = gdal.Open(item)  
            im_width = dataset.RasterXSize 
            im_height = dataset.RasterYSize  
            im_geotrans = dataset.GetGeoTransform()  
            im_proj = dataset.GetProjection()  
            im_data = dataset.ReadAsArray(0, 0, im_width, im_height) 
            image_array1 = np.array([image_array1, im_data])
            del dataset
            return im_geotrans, im_proj, image_array1
#数组均值缺失值替换计算，好像有个数组参数skipna跳过缺失值处理
def calculate_imagary_mean(data_array):
    data_array[data_array<0]=np.nan
    data_array = np.nanmean(data_array, axis=0)
    return data_array
# 写入影像
def write_img(filename, im_geotrans, im_proj, im_data):
        """
        pass
        """  
        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
        else:
            im_bands, (im_height, im_width) = 1, im_data.shape
       
        driver = gdal.GetDriverByName("GTiff")  
        dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

        dataset.SetGeoTransform(im_geotrans)  
        dataset.SetProjection(im_proj)  

        if im_bands == 1:
            dataset.GetRasterBand(1).WriteArray(im_data)  
        else:
            for i in range(im_bands):
                dataset.GetRasterBand(i + 1).WriteArray(im_data[i])

        del dataset  

if __name__ == "__main__":
    os.chdir(r'path')
    geotrans, proj, image_array = read_tif_file(r'path') 
    image_tif = calculate_imagary_mean(image_array)
    image_result = write_img (r'path\Result.tif', geotrans, proj, image_tif)