python遥感影像最大值合成

前言

由于数据处理的区域受云的影响较大，经过去云操作后出现大量0值，故想用该影像周围几天的像元值来填补有云影像

---

---

代码

一、引入库？

from osgeo import gdal
import numpy as np
import os

二、完整代码

1.函数

def write_tif(data_raw,targetdata,output_name):
    driver = gdal.GetDriverByName('GTiff')
    out_file = driver.Create(output_name,targetdata.shape[1],targetdata.shape[0],1,6)
    out_file.SetProjection(data_raw.GetProjection())
    out_file.SetGeoTransform(data_raw.GetGeoTransform())
    out_file.GetRasterBand(1).WriteArray(targetdata)
    del out_file

2.主体

input_file = 'G:/Test/'
output_file = 'G:/Test/Test_mvc/'
prefix = '.tif'

if not os.path.exists(output_file):
    os.mkdir(output_file)

file_all = os.listdir(input_file)
file_num = len(file_all)

for i in range(file_num):
    if file_all[i].endswith(prefix):
        data_box = ''
        interpol_data = gdal.Open(input_file + file_all[i])
        interpoled_data = interpol_data.ReadAsArray()
        output_name = output_file + os.path.splitext(file_all[i])[0] + '_mvc.tif'
        if i == 0:
            data_box = file_all[i:i+2]
            for j in range(len(data_box)):
                if file_all[i] == data_box[j]:
                    continue
                else:
                    data = gdal.Open(input_file + data_box[j])
                    data_array = data.ReadAsArray()
                    interpoled_data = (data_array >= interpoled_data)*data_array + (interpoled_data > data_array)*interpoled_data
            write_tif(data,interpoled_data,output_name)
        elif i == 1:
            data_box = file_all[i-1:i+2]
            for j in range(len(data_box)):
                if file_all[i] == data_box[j]:
                    continue
                else:
                    data = gdal.Open(input_file + data_box[j])
                    data_array = data.ReadAsArray()
                    interpoled_data = (data_array >= interpoled_data)*data_array + (interpoled_data > data_array)*interpoled_data
            write_tif(data,interpoled_data,output_name)
        elif i == file_num-1:
            data_box = file_all[i-2:i+1]
            for j in range(len(data_box)):
                if file_all[i] == data_box[j]:
                    continue
                else:
                    data = gdal.Open(input_file + data_box[j])
                    data_array = data.ReadAsArray()
                    interpoled_data = (data_array >= interpoled_data)*data_array + (interpoled_data > data_array)*interpoled_data
            write_tif(data,interpoled_data,output_name)
        elif i == file_num:
            data_box = file_all[i-2:i+1]
            for j in range(len(data_box)):
                if file_all[i] == data_box[j]:
                    continue
                else:
                    data = gdal.Open(input_file + data_box[j])
                    data_array = data.ReadAsArray()
                    interpoled_data = (data_array >= interpoled_data)*data_array + (interpoled_data > data_array)*interpoled_data
            write_tif(data,interpoled_data,output_name)
        else:
            data_box = file_all[i-2:i+2]
            for j in range(len(data_box)):
                if file_all[i] == data_box[j]:
                    continue
                else:
                    data = gdal.Open(input_file + data_box[j])
                    data_array = data.ReadAsArray()
                    interpoled_data = (data_array >= interpoled_data)*data_array + (interpoled_data > data_array)*interpoled_data
            write_tif(data,interpoled_data,output_name)

结果对比

此处以处理前后的同一天影像为例

总结

以上就是以近5天为窗口对整个文件夹的影像进行滑动最大值合成的代码。