其实问题的核心就只是给一个矩阵一个 EPSG 的信息
一般来说，深度学习之后的影像是没有投影坐标系的，如果这时候想用软件（ENVI、ArcGIS）去拼接就不现实了，后来想着本科的时候写过一个 MATLAB 的代码，就是植被覆盖度的逐像素回归，当时用 MATLAB 做的，也在那篇 blog 写过，其实处理遥感问题，拿出投影信息和仿射信息然后处理矩阵就好了。即，我有一个有有投影坐标系的 RGB 三通的图像，还有一个预测的 mask 图像，我只需要读进 mask 的矩阵，然后用 gdal
或者 geotiffread
读出 RGB 的投影和仿射信息，然后写给 mask 就可以了（注意的是，dir 下面的文件名必须严格对应，要不然，矩阵都飞走了，对应地图上的位置都错了，还怎么去正确拼接）。
废话不多说，先上代码...
# -*- coding：utf-8
# __author__:Neverland
# __date__:
from osgeo import gdal
from PIL import Image
import numpy as np
import os
# initialize
outband_size = 1
# 原始影像
image_dir = './Bing-input/'
image_files = os.listdir(image_dir)
image_files.sort(key= lambda x:int(x[:-4]))
# 预测图像
mask_dir = './Bing-predict-mask/'
# ref-mask
ref_mask_dir = './Bing-Ref/'
# 打开影像
for i in range(len(image_files)):
#设置文件名
image_file = os.path.join(image_dir,image_files[i])
mask_file = os.path.join(mask_dir,image_files[i])
file = os.path.join(ref_mask_dir,image_files[i])
# 打开 RGB 影像获取信息
in_ds = gdal.Open(image_file)
width = in_ds.RasterXSize
height = in_ds.RasterYSize
#投影信息和仿射变换信息
in_ds_projection = in_ds.GetProjection()
in_ds_Transform = in_ds.GetGeoTransform()
#读取一个波段即可，主要是为了获取它的 Datatype，注意要用 DataType，不要用 dtype，numpy.array 没有 DataType
in_band1 = in_ds.GetRasterBand(1)
in_band1_array = in_band1.ReadAsArray()
datatype = in_band1.DataType
# 读取 mask 的矩阵
mask_image = Image.open(mask_file)
mask_array = np.asarray(mask_image)
mask_image.close()
#创建驱动
gtif_driver = gdal.GetDriverByName("GTiff")
out_ds = gtif_driver.Create(file,width,height,outband_size,datatype)
out_ds.SetProjection(in_ds_projection)
out_ds.SetGeoTransform(in_ds_Transform)
out_ds.GetRasterBand(1).WriteArray(mask_array)
out_ds.FlushCache()
del out_ds,mask_array
思路很简单，就是用 gdal
读取 RGB 的投影和仿射信息，然后用 PIL
读取 mask 的像素矩阵，然后写出即可。
from osgeo import gdal
from PIL import Image
import numpy as np
import os
导入所需要的包，gdal
估计有点难装，用起来语法还很奇怪。
image_dir = './Bing-input/'
image_files = os.listdir(image_dir)
image_files.sort(key= lambda x:int(x[:-4]))
# 预测图像
mask_dir = './Bing-predict-mask/'
# ref-mask
ref_mask_dir = './Bing-Ref/'
这里设置了 RGB 的文件的目录，然后用 lambda
表达式对数字重新排序，我以前在做 yolo 处理的时候直接 pop
是乱序的，之前发现了这个问题，所以，特地用这个方法先排序，当然也可以用这个方法重新的批量重命名。
然后设置 mask 的 dir，设置输出图像的路径。
for i in range(len(image_files)):
loop code
这个大循环，主要是为了来控制我的文件名，然后打开 RGB、打开 mask，写出的。
image_file = os.path.join(image_dir,image_files[i])
mask_file = os.path.join(mask_dir,image_files[i])
file = os.path.join(ref_mask_dir,image_files[i])
这里就是用 os.path.join(dir,file_name)
初始化 RGB、mask 还有 写出文件的文件名的。
# 打开 RGB 影像获取信息
in_ds = gdal.Open(image_file)
width = in_ds.RasterXSize
height = in_ds.RasterYSize
#投影信息和仿射变换信息
in_ds_projection = in_ds.GetProjection()
in_ds_Transform = in_ds.GetGeoTransform()
#读取一个波段即可，主要是为了获取它的 Datatype，注意要用 DataType，不要用 dtype，numpy.array 没有 DataType
in_band1 = in_ds.GetRasterBand(1)
in_band1_array = in_band1.ReadAsArray()
datatype = in_band1.DataType
首先 gdal.Open
打开 RGB 文件，获取长和宽，然后得到投影和仿射信息，然后读取一个波段即可，为的是得到它的 DataType
，不要以为它是个 numpy.array
就尝试用 dtype
去获取它的数据类型，gdal
有专门的数据类型，如果 numpy.dtype
只是 uint8
而已。
mask_image = Image.open(mask_file)
mask_array = np.asarray(mask_image)
mask_image.close()
这里用 Image.io
读取 mask，然后转成 numpy
的矩阵，这时候文件就可以关掉，没必要占内存了。
gtif_driver = gdal.GetDriverByName("GTiff")
out_ds = gtif_driver.Create(file,width,height,outband_size,datatype)
out_ds.SetProjection(in_ds_projection)
out_ds.SetGeoTransform(in_ds_Transform)
out_ds.GetRasterBand(1).WriteArray(mask_array)
out_ds.FlushCache()
del out_ds,mask_array
创建个 GTiff
的驱动，并且驱动创建要写的文件的信息（很不喜欢这种格式，一堆记不住的代码，经常要回去复制），然后 SetProjection
和 SetGeoTransform
然后再写出矩阵，磁盘缓存，最后，记得销毁 out_ds
和 mask_array
的内存，保不准你的矩阵就太大卡巴斯基了
呢?
Python 代码就到这里了...有点长还有点啰嗦，不大好记。
# -*- coding：utf-8
# __author__:Neverland
# __date__:
image_dir = './Bing-input/';
mask_dir = './Bing-predict-mask/';
output_dir = './Bing-Ref/';
for i = 1:63
%拼接文件名
image_name = strcat(num2str(i),'.tif');
image_path = strcat(image_dir,image_name)
mask_path = strcat(mask_dir,image_name)
%读取矩阵，写矩阵
[A,R] = geotiffread(image_path);
info = geotiffinfo(image_path);
mask_image = imread(mask_path);
output_name = strcat(output_dir,image_name);
geotiffwrite(output_name,mask_image,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
end
相比之下，MATLAB 真就简单太多了，除了最后一句代码记不住那几个结构体，其他简单通俗易懂，最后一句直接去以前的 CSDN 复制一下，改下矩阵的名称，和输出名的变量名就好了，就不详细解释了。
拼接的代码我倒还没去写，直接 ENVI 拼然后导出的，但是估计尝试个 1 个小时应该也能写出来，一般写代码就是从少到多，比如拼两张，再到一个 for
循环全部拼起来。
不要问我为什么那么久没更了，天天跟我导打电话就好久了。。
性情中人点首歌：
