Python遥感开发之批量拼接
- 1 遥感图像无交错的批量拼接
- 2 遥感图像有交错的批量拼接
前言:主要借助python实现遥感影像的批量拼接,遥感影像的批量拼接主要分为两种情况,一种是遥感图像无交错,另一种情况是遥感图像相互有交错。具体实现请参考以下代码,如有问题请及时反馈。
1 遥感图像无交错的批量拼接
此方法是各个遥感文件是没有相互交错的拼接,如下图所示。个人可以使用Arcgis进行查看。
实现思路:通过每个遥感数据的经纬度进行拼接下一个遥感数据文件。
import os
from osgeo import gdal
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 + "\\" + f
out_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 list1
def 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_list02
def 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.RasterXSize
rows = ds.RasterYSize
ingeo = ds.GetGeoTransform()
proj = ds.GetProjection()
minx = ingeo[0]
maxy = ingeo[3]
maxx = ingeo[0] + ingeo[1] * cols
miny = ingeo[3] + ingeo[5] * rows
ds = None
for file in infile_list[1:]:
ds = gdal.Open(file)
cols = ds.RasterXSize
rows = ds.RasterYSize
geo = ds.GetGeoTransform()
minx_ = geo[0]
maxy_ = geo[3]
maxx_ = geo[0] + geo[1] * cols
miny_ = geo[3] + geo[5] * rows
minx = min(minx, minx_)
maxy = max(maxy, maxy_)
maxx = max(maxx, maxx_)
miny = min(miny, miny_)
geo = None
ds = None
newcols = 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 = None
outband.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 os
def 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 + "\\" + f
out_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 list1
def 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 data
def get_same_list(image, infile_list):
infile_list02 = []
for data in infile_list:
if image in data:
infile_list02.append(data)
return infile_list02
if __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("-------")