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