Python遥感开发之批量拼接

前言：主要借助python实现遥感影像的批量拼接，遥感影像的批量拼接主要分为两种情况，一种是遥感图像无交错，另一种情况是遥感图像相互有交错。具体实现请参考以下代码，如有问题请及时反馈。

---

1 遥感图像无交错的批量拼接

此方法是各个遥感文件是没有相互交错的拼接，如下图所示。个人可以使用Arcgis进行查看。

实现思路：通过每个遥感数据的经纬度进行拼接下一个遥感数据文件。

import os
from osgeo import gdaldef get_data_list(file_path, out = ""):list1 = []  # 文件的完整路径if os.path.isdir(file_path):fileList = os.listdir(file_path)if out != "":for f in fileList:out_data = out + "\\" + fout_data = out_data.replace(".HDF", "_ndvi.tif")list1.append(out_data)else:for f in fileList:pre_data = file_path + '\\' + f  # 文件的完整路径list1.append(pre_data)return list1def get_same_list(image, infile_list):infile_list02 = []for data in infile_list:if image in data:# print("----", data)infile_list02.append(data)return infile_list02def get_same_image_list(infile_list):image_list= []for file in infile_list:filename = file[-31:-23]if filename not in image_list:image_list.append(filename)return list(set(image_list))def pinjie(infile_list,outfile):ds = gdal.Open(infile_list[0])cols = ds.RasterXSizerows = ds.RasterYSizeingeo = ds.GetGeoTransform()proj = ds.GetProjection()minx = ingeo[0]maxy = ingeo[3]maxx = ingeo[0] + ingeo[1] * colsminy = ingeo[3] + ingeo[5] * rowsds = Nonefor file in infile_list[1:]:ds = gdal.Open(file)cols = ds.RasterXSizerows = ds.RasterYSizegeo = ds.GetGeoTransform()minx_ = geo[0]maxy_ = geo[3]maxx_ = geo[0] + geo[1] * colsminy_ = geo[3] + geo[5] * rowsminx = min(minx, minx_)maxy = max(maxy, maxy_)maxx = max(maxx, maxx_)miny = min(miny, miny_)geo = Noneds = Nonenewcols = int((maxx - minx) / abs(ingeo[1]))newrows = int((maxy - miny) / abs(ingeo[5]))driver = gdal.GetDriverByName("GTiff")outds = driver.Create(outfile, newcols, newrows, 1, gdal.GDT_Int16)outgeo = (minx, ingeo[1], 0, maxy, 0, ingeo[5])outds.SetGeoTransform(outgeo)outds.SetProjection(proj)outband = outds.GetRasterBand(1)for file in infile_list:ds = gdal.Open(file)data = ds.ReadAsArray()geo = ds.GetGeoTransform()x = int(abs((geo[0] - minx) / ingeo[1]))y = int(abs((geo[3] - maxy) / ingeo[5]))outband.WriteArray(data, x, y)ds = Noneoutband.FlushCache()pass
if __name__ == '__main__':infile = r"C:\Users\Administrator\Desktop\01提取ndvi"outfile = r"C:\Users\Administrator\Desktop\02拼接"infile_list = get_data_list(infile)image_name_list = get_same_image_list(infile_list)print(image_name_list)for name in image_name_list:print(name)infile_list02 = get_same_list(name, infile_list)pinjie(infile_list02,outfile+"\\"+name+".tif")

2 遥感图像有交错的批量拼接

此方法是各个遥感文件是有相互交错的拼接，如下图所示，具体可以使用Arcgis进行查看。

实现思路：借助gdal中WarpOptions的方法实现，有点类似于镶嵌

import numpy as np
from osgeo import gdal, gdalconst
import osdef RasterMosaic(firstinputfilePath, inputfileList, outputfilePath):inputrasfile1 = gdal.Open(firstinputfilePath, gdal.GA_ReadOnly)  # 第一幅影像inputProj1 = inputrasfile1.GetProjection()options = gdal.WarpOptions(srcSRS=inputProj1, dstSRS=inputProj1, format='GTiff')gdal.Warp(outputfilePath, inputfileList, options=options)def get_data_list(file_path, out=""):list1 = []  # 文件的完整路径if os.path.isdir(file_path):fileList = os.listdir(file_path)if out != "":for f in fileList:out_data = out + "\\" + fout_data = out_data.replace(".HDF", "_ndvi.tif")list1.append(out_data)else:for f in fileList:pre_data = file_path + '\\' + f  # 文件的完整路径list1.append(pre_data)return list1def get_same_image_list(infile_list):image_list = []for file in infile_list:filename = file[-20:-12]if filename not in image_list:image_list.append(filename)return list(set(image_list))def get_infile(image,infile_list):for data in infile_list:if image in data:return datadef get_same_list(image, infile_list):infile_list02 = []for data in infile_list:if image in data:infile_list02.append(data)return infile_list02if __name__ == '__main__':inputfile_path = r"D:\风云数据\MERSI-II陆表反射比1KM段产品\b1\01原始"outfile = r"D:\风云数据\MERSI-II陆表反射比1KM段产品\b1\02拼接"infile_list = get_data_list(inputfile_path)image_list = get_same_image_list(infile_list)print(image_list)for image in image_list:firstinputfilePath = get_infile(image,infile_list)infile_list02 = get_same_list(image, infile_list)print(image)print(firstinputfilePath)print(infile_list02)RasterMosaic(firstinputfilePath, infile_list02, outfile+"\\"+image+"_b1.tif")print("-------")