核心就是找到位置赋值，然后写仿射信息和投影信息！
Export Layer to TIFF 会导出三个波段，在 ArcGIS 里显示很不方便，要用 ENVI 的
Save as功能保存成单波段的图片，但是吧，软件用起来保不齐哪天因为矩阵太大就
卡巴斯基了，所以在网上找了一圈代码，发现用不了，看了眼，自己想了想思路，把思路理清了，就自己写了个代码，虽然可能有重复的执行代码部分，但是也就是自己用着玩玩而已，毕竟一句
python combine.py就把一个文件夹里面的文件都拼接了。能够极大提升自己代码效率，还有个好处，我个人自己认为的:
第一种方法导出图层，导出了3个波段图像，第二种方法保存图像，保存成一个单波段的图像，还是可以的，但是说实话。
打开好多图然后慢慢生成 *.enp
再用 Seamless Mosaic
不知道要多久，代码处理这种小范围的可能 3 秒钟就解决战斗了...
from osgeo import gdal
import os
import numpy as np
import matplotlib.pyplot as plt
"""
# input-arg:输入图像的文件夹
获取所有图像的尺寸，并得到需要创建矩阵的尺寸
"""
def get_size(image_dir):
#
#排序
image_files = os.listdir(image_dir)
image_files.sort(key = lambda x:int(x[:-4]))
#初始化
row = 0
col = 0
init_img = gdal.Open(os.path.join(image_dir,image_files[0]))
#第一块图像的行列
row += init_img.RasterXSize
col += init_img.RasterYSize
#获取初始值的 left_X,left_y，即左上角的点
init_img_transform = init_img.GetGeoTransform()
init_left_x = init_img_transform[0]
init_left_y = init_img_transform[3]
width_pixel = init_img_transform[1]
height_pixel = abs(init_img_transform[5])
del init_img
for file in image_files[1:]:
#就直接拼接文件名，丑就丑吧，自己写的代码还想咋样呢...
in_ds = gdal.Open(os.path.join(image_dir,file))
#获取仿射矩阵的信息，然后和第一块去比较判断
in_ds.transform = in_ds.GetGeoTransform()
#仿射矩阵的左上角的 x,y
left_x = in_ds.transform[0]
left_y = in_ds.transform[3]
#在第一行或者第一列时才计数
if (left_x != init_left_x) and (left_y == init_left_y):
col += in_ds.RasterXSize
elif (left_y != init_left_y) and (left_x == init_left_x):
row += in_ds.RasterYSize
del in_ds
return row,col
def give_data(image_dir,temp_array):
image_files = os.listdir(image_dir)
image_files.sort(key = lambda x:int(x[:-4]))
#初始化
row = 0
col = 0
init_img = gdal.Open(os.path.join(image_dir,image_files[0]))
#获取关键信息
init_img_transform = init_img.GetGeoTransform()
init_left_x = init_img_transform[0]
init_left_y = init_img_transform[3]
width_pixel = init_img_transform[1]
height_pixel = abs(init_img_transform[5])
#初始化的赋值
temp_array[0:init_img.RasterYSize,0:init_img.RasterXSize] = init_img.GetRasterBand(1).ReadAsArray()
for file in image_files[1:]:
in_ds = gdal.Open(os.path.join(image_dir,file))
in_ds.transform = in_ds.GetGeoTransform()
left_x = in_ds.transform[0]
left_y = in_ds.transform[3]
#计算他们所在的行列，但是由于裁剪的时候是向下裁剪，所以要用第一张图片的 y 减去后面的图片
new_row = int((init_left_y - left_y)/(height_pixel))
new_col = int((left_x - init_left_x)/(width_pixel))
#找到在矩阵索引的位置
img_x_loc = new_row + in_ds.RasterYSize
img_y_loc = new_col + in_ds.RasterXSize
temp_array[new_row:img_x_loc,new_col:img_y_loc] = in_ds.GetRasterBand(1).ReadAsArray()
del in_ds
return temp_array
#写出tif
def write_array(image_dir,temp_array):
image_files = os.listdir(image_dir)
init_img = gdal.Open(os.path.join(image_dir,image_files[0]))
img_Projection = init_img.GetProjection()
img_GeoTransform = init_img.GetGeoTransform()
DataType = init_img.GetRasterBand(1).DataType
width = temp_array.shape[1]
height = temp_array.shape[0]
gtif_driver = gdal.GetDriverByName("GTiff")
out_file_name = './combine.tif'
outband_size = 1
out_ds = gtif_driver.Create(out_file_name,width,height,outband_size,DataType)
#投影和仿射信息直接用第一幅图像的信息即可
out_ds.SetProjection(img_Projection)
out_ds.SetGeoTransform(img_GeoTransform)
out_ds.GetRasterBand(1).WriteArray(temp_array)
out_ds.FlushCache()
del out_ds
# main 测试
if __name__ == '__main__':
mask_dir = './Bing-Ref/'
#获取文件夹中的图片，计算他们拼接需要的尺寸和内存
[row,col] = get_size(mask_dir)
#创建一个临时的虚拟矩阵，用来接收图片的像素
temp = np.zeros([row,col])
#一个一个打开文件然后一个一个赋值到矩阵的相应位置
temp = give_data(mask_dir,temp)
#写出带参考系的矩阵
write_array(mask_dir,temp)
先上代码，同时，这个是我自己写的，所以会有些冗余代码，而且前期的要求挺高的。
1.要用我以前说的那种裁剪方法，因为这个拼接方法是那种裁剪的反序。2.文件夹中没有冗余文件，要不文件名排序的时候挂了我也没办法（自己懒，不想写太多的错误判断）3.文件名自己命名，1,2.3.4.5.6...讲道理其实矩阵不连续也没问题，因为有根据仿射矩阵判断位置的功能。
至于文件名怎么批量重新命名，以后有空再写篇文章讲一讲，要不手动删或者对应比较废时间。
不想用 PPT 画图了，用 ENVI 的拼接功能做个演示，偷偷写篇博客，然后赶 DDL 去了，但是这篇估计得写 1 小时左右。
如上图所示，用 gdal.Open('file')
之后再使用 file.GetGeoTransform()
之后会得到一个仿射矩阵的信息，其中第 1 个元素（0起始）代表的是左上角的 x 坐标
,第 4 个元素是 左上角的 y 坐标
，然后第 2 个元素是东西向分辨率，第 6 个元素是 南北向分辨率（负值，要取绝对值）。
代码的思路也很简单，就是通过计算两幅图像之间的 X,Y 坐标差，然后除以他们各自方向上的分辨率，计算他们在合并矩阵的起始位置，然后再合并矩阵的位置的 [start_row:start_row + img.RasterYSzie,start_col:start_col + img.RasterXSize] 赋上相应的图像的值即可。
在细分下代码的思路，分为 3步，3个函数：
1.得到合并矩阵的 shape
，即 get_size
函数，在文件夹中自动计算。2.根据 left_x
,left_y
自动计算合并矩阵的位置，然后赋值。3.由于合并矩阵的 geotransform
和第一张图一样，所以我们只要把 投影信息 和仿射矩阵赋给合并的矩阵即可。
def get_size(image_dir):
#
#排序
image_files = os.listdir(image_dir)
image_files.sort(key = lambda x:int(x[:-4]))
#初始化
row = 0
col = 0
init_img = gdal.Open(os.path.join(image_dir,image_files[0]))
#第一块图像的行列
row += init_img.RasterXSize
col += init_img.RasterYSize
#获取初始值的 left_X,left_y，即左上角的点
init_img_transform = init_img.GetGeoTransform()
init_left_x = init_img_transform[0]
init_left_y = init_img_transform[3]
width_pixel = init_img_transform[1]
height_pixel = abs(init_img_transform[5])
del init_img
for file in image_files[1:]:
#就直接拼接文件名，丑就丑吧，自己写的代码还想咋样呢...
in_ds = gdal.Open(os.path.join(image_dir,file))
#获取仿射矩阵的信息，然后和第一块去比较判断
in_ds.transform = in_ds.GetGeoTransform()
#仿射矩阵的左上角的 x,y
left_x = in_ds.transform[0]
left_y = in_ds.transform[3]
#在第一行或者第一列时才计数
if (left_x != init_left_x) and (left_y == init_left_y):
col += in_ds.RasterXSize
elif (left_y != init_left_y) and (left_x == init_left_x):
row += in_ds.RasterYSize
del in_ds
return row,col
get_size函数吧，输入是
mask_dir，然后首先读取左上角的图像做一个初始化，然后得到相应的初始值，得到信息之后记得关闭内存，别占着呗。然后一个一个打开后续的文件，然后得到仿射矩阵的
left_x,
left_y然后最重要的一件事是下面的条件的判断，即，我所在的图像在 第一行 行或者第一列才要加到相应的行列数，要不然你会得到一个很奇特的数字，
col对应的是列，是图像的宽，对应
RasterXSize，
row对应的是行，对应的是图像的高度，对应
RasterYSize。然后，这样就得到合并矩阵的尺寸了，然后返回
行和
列数，在主函数的地方可以用
numpy去创一个临时的矩阵，然后用来接各个图像的值了。
def give_data(image_dir,temp_array):
image_files = os.listdir(image_dir)
image_files.sort(key = lambda x:int(x[:-4]))
#初始化
row = 0
col = 0
init_img = gdal.Open(os.path.join(image_dir,image_files[0]))
#获取关键信息
init_img_transform = init_img.GetGeoTransform()
init_left_x = init_img_transform[0]
init_left_y = init_img_transform[3]
width_pixel = init_img_transform[1]
height_pixel = abs(init_img_transform[5])
#初始化的赋值
temp_array[0:init_img.RasterYSize,0:init_img.RasterXSize] = init_img.GetRasterBand(1).ReadAsArray()
for file in image_files[1:]:
in_ds = gdal.Open(os.path.join(image_dir,file))
in_ds.transform = in_ds.GetGeoTransform()
left_x = in_ds.transform[0]
left_y = in_ds.transform[3]
#计算他们所在的行列，但是由于裁剪的时候是向下裁剪，所以要用第一张图片的 y 减去后面的图片
new_row = int((init_left_y - left_y)/(height_pixel))
new_col = int((left_x - init_left_x)/(width_pixel))
#找到在矩阵索引的位置
img_x_loc = new_row + in_ds.RasterYSize
img_y_loc = new_col + in_ds.RasterXSize
temp_array[new_row:img_x_loc,new_col:img_y_loc] = in_ds.GetRasterBand(1).ReadAsArray()
del in_ds
return temp_array
第二步就是给矩阵相应位置附上所需要的的值了，give_data
输入参数是图像的文件夹，然后输出一个赋值好的矩阵。为什么说代码有重复呢，就是因为第一步和第二步很多都是重复的，但是我只是比较懒，不想去多修改了。 具体思路就是之前所说的，首先根据图像的 left_x
和 left_y
（减或者反减) init_left_x
和 init_left_y
，然后除以他们各自方向的分辨率，得到他们的起始位置的 row,col
，再用 row+img.RasterYSize
,和 col+img.RasterXSize
计算所在的行列数，并对相应的位置赋值
plus:吐槽下，一会儿用 MATLAB 一会儿用 Python 真的很难受，矩阵的索引老是习惯的 ()
和 []
分不清。
最后返回赋值完成的矩阵。
def write_array(image_dir,temp_array):
image_files = os.listdir(image_dir)
init_img = gdal.Open(os.path.join(image_dir,image_files[0]))
img_Projection = init_img.GetProjection()
img_GeoTransform = init_img.GetGeoTransform()
DataType = init_img.GetRasterBand(1).DataType
width = temp_array.shape[1]
height = temp_array.shape[0]
gtif_driver = gdal.GetDriverByName("GTiff")
out_file_name = './combine.tif'
outband_size = 1
out_ds = gtif_driver.Create(out_file_name,width,height,outband_size,DataType)
#投影和仿射信息直接用第一幅图像的信息即可
out_ds.SetProjection(img_Projection)
out_ds.SetGeoTransform(img_GeoTransform)
out_ds.GetRasterBand(1).WriteArray(temp_array)
out_ds.FlushCache()
del out_ds
plus:如果真的老是分不清 col
、row
、img.RasterXSize
和 img.RasterYSize
还有 array.shape[0]
和 array.shape[1]
的话要不做个表要不就试试吧。
然后检查下，图像和 ENVI 拼接保存的图像，发现文件大小竟然一样，看来思路可能是一样的呢。。主要是看看有没有因为上面这个表格老是记不清位置写错代码导致赋值错地方或者导致矩阵的尺寸不对。
| 宽 | 高 |
| row | col |
| img.RasterYSize | img.RasterXSize |
| array.shape[1] | array.shape[0] |
文章转载自狗头Rser，如果涉嫌侵权，请发送邮件至：contact@modb.pro进行举报，并提供相关证据，一经查实，墨天轮将立刻删除相关内容。
