Python遥感影像拼接

对于栅格影像，我们一般可以采用ENVI或ArcGIS平台进行拼接，也可以通过GEE或PIE Engine云平台进行处理。如果我们想利用人工神经网络进行操作，就需要云平台中的数据导出到本地。对于大尺度的遥感影像而言，就会存在每一景影像大小受限的问题，就会分成许多景影像，如下图所示：

因此，在对其进行操作之前，我们需要对其进行拼接。接下来我们就介绍如何利用Python进行栅格影像的拼接。

在Python中有两个强大的模块，一个是raster，一个是gdal，二者都可以对栅格数据进行处理：

Raster

import rasterio
from rasterio.merge import merge
from rasterio.plot import show
import glob
import os
import matplotlib.pyplot as plt

Importing required modules and find all `tif` files from the folder

# File and folder paths
dirpath = r"E:\数据集\祁连山L8"
out_fp = os.path.join(dirpath, "qilianshanL8.tif")
tif_file = glob.glob(os.path.join(dirpath, "*.tif"))
print(tif_file)

create a list for the source raster datafiles

src_files_to_mosaic = []
for tif_f in tif_file:src = rasterio.open(tif_f)src_files_to_mosaic.append(src)
print('src_files_to_mosaic', src_files_to_mosaic)

update the metadata with our new dimensions, transform and CRS and write to our computer

out_meta.update({"driver": "GTiff","height": mosaic.shape[1],"width": mosaic.shape[2],"transform": out_trans,"crs": "EPSG:4326"})
# Write the mosaic raster to disk
with rasterio.open(out_fp, "w", **out_meta) as dest:dest.write(mosaic)

在处理完之后，我们也可以plt出拼接影像的效果图

mosaic, out_trans = merge(src_files_to_mosaic)
# Plot the result
show(mosaic, cmap='terrain')

当然，需要注意的是，在拼接过程中，参与拼接的影像的投影参数应该保持一致，我们既可以手动输入，也可以采用自动计算的方式

"crs": "+proj=utm +zone=35 +ellps=GRS80 +units=m +no_defs "
"crs": "EPSG:4326"

但是，对于大尺度遥感影像而言，分幅影像的投影参数往往是不同的，因此我们先要将其投影参数统一化，然后再进行拼接操作。

for i in range(len(tif_file)):dstfilename = "EPGSG32649" + str(i) + ".tif" # 根据自己的需求设置文件名with rasterio.open(tif_file[i]) as src:transform, width, height = calculate_default_transform(src.crs, dst_crs, src.width, src.height, *src.bounds)kwargs = src.meta.copy()kwargs.update({"crs": dst_crs, "transform": transform, "width": width, "height": height, "compress":'lzw'})with rasterio.open(dstfilename, "w", **kwargs) as dst:for i in range(1, src.count + 1):reproject(source=rasterio.band(src, i),destination=rasterio.band(dst, i),src_transform=src.transform,src_crs=src.crs,dst_transform=transform,dst_crs=dst_crs,resampling=Resampling.nearest,)

下面是官方关于Creating a raster mosaic的链接：

Rasterhttps://autogis-site.readthedocs.io/en/latest/notebooks/Raster/raster-mosaic.html

GDALGDALhttps://www.neonscience.org/resources/learning-hub/tutorials/merge-lidar-geotiff-py