在使用 GDAL 读取遥感影像的 dateset 后，通过 Projection 和 GeoTransform 这两个参数即可确定影像的地理位置。对于开源的遥感影像数据集，它们的这两个主要参数一般都会进行“脱敏”处理，即去除投影坐标系信息，修饰 GeoTransform 的范围，从而使影像地理坐标位于笛卡尔坐标系下。去除投影信息的结果是显然意见的，那么，对六个 GeoTransform 参数进行修饰会带来怎样的影响呢？一次批量裁剪数据时出现的影响错位问题困扰了我许久，拨开云雾查明原因后，才发现这个问题与这些参数密切相关......

简介 GeoTransform

在遥感图像中，x 与 y 分别对应行与列，原点对应于左上角的像素坐标。使用 gdal 读取 GeoTransform 的方式如下：

from osgeo import gdal
filename = "D:\\Desktop\a.tif"
ds = gdal.Open(filename)
geo_trans = ds.GetGeoTransform()
print(geo_trans)
del ds

对于未做特殊处理的遥感影像，结果一般如下所示：

(486892.5, 15.0, 0.0, 4105507.5, 0.0, -15.0)

这六个参数的详细解释如下图所示：

其中需留意的是，这里的 geo_trans[5] 为负数，也正因为此，图像才能从左上角的 (0, 0) 铺展到右下角。如果这个值脱离一般情形，变为正数的话，那么结果会怎样呢？

反面教材

以下为反面实例。如果不调整 geo_trans[5] 的话，裁剪出来的影像将如下图，虽然保留了正确的图像信息，但位置发生了偏移。

打开 ArcGIS 查看原始图像的基本属性：

从中可以看到，坐标系已做了脱敏处理。试着用 GDAL 来读取一下 geotrans 试试：

发现 geo_trans 也进行了脱敏处理。留意这里显示的左上角坐标与 ArcGIS 中显示的有所不同。就是这两个小细节会影响最终结果的整体表现，所以处理数据一定要细心！

改进后的随机裁剪

随机裁剪常用于批量生成训练集，测试算法性能等。一副遥感影像可以被 GDAL 整体读取为 np.array，随机裁剪正可以由此进行。一个简单的思路是从左上方开始，引入 random 库随机在可行的位置定位小窗口的左上角，然后进行裁剪，按坐标偏移量更新 geotrans，最后对裁剪出的窗口设置与原图像相同的投影并导出。

在上述的基础上，我额外进行了如下改进：

• 引入全局变量记录左上角坐标值以避免裁剪出重复图像；
• 对脱敏数据额外进行判断，避免裁剪出的图像出现坐标偏移；
• 结合实际数据加入 +—0.5 的坐标偏离量（非必需，具体问题具体分析）；
• 确保裁剪出图像左上角坐标位置为 20 的整数倍，以匹配后续图像。

核心代码如下，已经写了简要的注释：

from osgeo import gdal
import numpy as np
import random


coords_done = []  # check duplicated situation

def random_crop(ImagePath, ImageSavePath, CropSize, CutNum):
    """
    ImagePath: The input raster to be clipped
    ImageSavePath: The path to save the clipped images
    CropSize: The size of each clipped image (unit: pixel)
    CutNum: The number of exported images with `CropSize`
    """
    dataset_img = readTif(ImagePath)
    width = dataset_img.RasterXSize
    height = dataset_img.RasterYSize
    proj = dataset_img.GetProjection()
    geotrans = dataset_img.GetGeoTransform()
    img = dataset_img.ReadAsArray(0, 0, width, height)
    
    geotrans_list = []
    geotrans_list.append(geotrans[0])
    geotrans_list.append(geotrans[1])
    geotrans_list.append(geotrans[2])
    geotrans_list.append(geotrans[3])
    geotrans_list.append(geotrans[4])
    if geotrans[5] < 0:
        geotrans_list.append(geotrans[5]*-1)
    else:
        geotrans_list.append(geotrans[5])

    if not os.path.exists(ImageSavePath):
        os.makedirs(ImageSavePath)

    fileNum = len(os.listdir(ImageSavePath))
    new_name = fileNum + 1
    
    global coords_done
    
    while (new_name < CutNum + fileNum + 1):
        # ensure the left upper coordinate is k*20 (k is integeral)
        UpperLeftX = 20 * random.randint(0, (height - CropSize) // 20)
        UpperLeftY = 20 * random.randint(0, (width - CropSize) // 20)
        if (UpperLeftX, UpperLeftY) not in coords_done:
            coords_done.append((UpperLeftX, UpperLeftY))
        else:
            print("Image Duplicate! Skip to next clipped image...")
            continue
        # When there is only one band, then `img.shape` will be 2 rather than 3
        if (len(img.shape) == 2):
            imgCrop = img[UpperLeftX: UpperLeftX + CropSize,
                      UpperLeftY: UpperLeftY + CropSize]
        else:
            imgCrop = img[:,
                      UpperLeftX: UpperLeftX + CropSize,
                      UpperLeftY: UpperLeftY + CropSize]
        x = int(UpperLeftX)  # row
        y = int(UpperLeftY)  # col
        # calculate new position
        px = geotrans[0] + y * geotrans[1] + x * geotrans[2]
        if geotrans[5] < 0:
            py = geotrans[3] + y * geotrans[4] + x * geotrans[5] * -1
        else:
            py = geotrans[3] + y * geotrans[4] + x * geotrans[5]
        # The +- 0.5 is needed to match the image accurately
        geotrans_list[0] = px - 0.5
        geotrans_list[3] = py + 0.5
        geotrans_tuple = tuple(geotrans_list)
        writeTiff(imgCrop, geotrans_tuple, proj, ImageSavePath + "/%d.tif" % new_name)
        new_name = new_name + 1
    coords_done = []

最终取得了心仪的结果（原始底图调整了波段组合与透明度以强调对比）：