晉中網(wǎng)站建設(shè)公司百度我的訂單
Python遙感開發(fā)之批量拼接
- 1 遙感圖像無交錯(cuò)的批量拼接
- 2 遙感圖像有交錯(cuò)的批量拼接
前言:主要借助python實(shí)現(xiàn)遙感影像的批量拼接,遙感影像的批量拼接主要分為兩種情況,一種是遙感圖像無交錯(cuò),另一種情況是遙感圖像相互有交錯(cuò)。具體實(shí)現(xiàn)請(qǐng)參考以下代碼,如有問題請(qǐng)及時(shí)反饋。
1 遙感圖像無交錯(cuò)的批量拼接
此方法是各個(gè)遙感文件是沒有相互交錯(cuò)的拼接,如下圖所示。個(gè)人可以使用Arcgis進(jìn)行查看。
實(shí)現(xiàn)思路:通過每個(gè)遙感數(shù)據(jù)的經(jīng)緯度進(jìn)行拼接下一個(gè)遙感數(shù)據(jù)文件。
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 遙感圖像有交錯(cuò)的批量拼接
此方法是各個(gè)遙感文件是有相互交錯(cuò)的拼接,如下圖所示,具體可以使用Arcgis進(jìn)行查看。
實(shí)現(xiàn)思路:借助gdal中WarpOptions的方法實(shí)現(xiàn),有點(diǎn)類似于鑲嵌
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:\風(fēng)云數(shù)據(jù)\MERSI-II陸表反射比1KM段產(chǎn)品\b1\01原始"outfile = r"D:\風(fēng)云數(shù)據(jù)\MERSI-II陸表反射比1KM段產(chǎn)品\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("-------")