基于GEE生成的分幅遥感影像的镶嵌处理。
不足之处是：镶嵌之后是由多波段影像变成单波段影像。。。
在该链接代码”基于python实现.tif格式遥感影像的镶嵌处理（合并）_蓝色风沙的技术博客_51CTO博客“基础上实现批量处理。
代码如下：
from osgeo import gdal
import os
import glob
import math
def GetExtent(in_fn):
ds=gdal.Open(in_fn)
geotrans=list(ds.GetGeoTransform())
xsize=ds.RasterXSize
ysize=ds.RasterYSize
min_x=geotrans[0]
max_y=geotrans[3]
max_x=geotrans[0]+xsize*geotrans[1]
min_y=geotrans[3]+ysize*geotrans[5]
ds=None
return min_x,max_y,max_x,min_y
def RasterMosaic(in_files,outpath):
in_fn=in_files[0]
#获取待镶嵌栅格的最大最小的坐标值
min_x,max_y,max_x,min_y=GetExtent(in_fn)
for in_fn in in_files[1:]:
minx,maxy,maxx,miny=GetExtent(in_fn)
min_x=min(min_x,minx)
min_y=min(min_y,miny)
max_x=max(max_x,maxx)
max_y=max(max_y,maxy)
#计算镶嵌后影像的行列号
in_ds=gdal.Open(in_files[0])
geotrans=list(in_ds.GetGeoTransform())
width=geotrans[1]
height=geotrans[5]
columns=math.ceil((max_x-min_x)/width)
rows=math.ceil((max_y-min_y)/(-height))
in_band=in_ds.GetRasterBand(1)
driver=gdal.GetDriverByName('GTiff')
out_ds=driver.Create(outpath,columns,rows,1,in_band.DataType)
out_ds.SetProjection(in_ds.GetProjection())
geotrans[0]=min_x
geotrans[3]=max_y
out_ds.SetGeoTransform(geotrans)
out_band=out_ds.GetRasterBand(1)
#定义仿射逆变换
inv_geotrans=gdal.InvGeoTransform(geotrans)
#开始逐渐写入
for in_fn in in_files:
in_ds=gdal.Open(in_fn)
in_gt=in_ds.GetGeoTransform()
#仿射逆变换
offset=gdal.ApplyGeoTransform(inv_geotrans,in_gt[0],in_gt[3])
x,y=map(float,offset)
print(x,y)
trans=gdal.Transformer(in_ds,out_ds,[])#in_ds是源栅格，out_ds是目标栅格
success,xyz=trans.TransformPoint(False,0,0)#计算in_ds中左上角像元对应out_ds中的行列号
x,y,z=map(float,xyz)
print(x,y,z)
data=in_ds.GetRasterBand(1).ReadAsArray()
out_band.WriteArray(data,x,y)#x，y是开始写入时左上角像元行列号
del in_ds,out_band,out_ds
inputfile_path= r'D:\cities\input'
list_files=os.listdir(inputfile_path)
#print(list_files)
#批量后文件的命名及保存
targetFile = []
for i in list_files:
targetFile.append(i[:-26] + i[-4:])
targetFile=list(set(targetFile))
print(targetFile)
for j in range(len(targetFile)):
mergename = targetFile[j]
name_str = mergename[:-4]
if os.path.exists(mergename):
os.remove(mergename)
# print (name_str)
inputfileList=[]
for readPath in list_files:
if readPath.startswith(name_str):
inputfilepath = inputfile_path + "\\" + readPath
# print(inputfilepath)
inputfileList.append(inputfilepath)
print(inputfileList)
outpath= r'D:\cities\Inputimages\\'
RasterMosaic(inputfileList,outpath+mergename)
print("==========================>"+mergename+"完成<============================")