用gdal分割tif文件，保持分割后的小文件位置不变

遥感tif文件太大了（30G），想将他切分为多个小的分别进行处理，考虑道之后想要在重组的，但是看了一圈都没有说明切分后的位置问题，就边查边写了一个。先看看效果图


可以看到切分后再将所有小块tif文件导入Arcgis位置互相没有偏差，黑框应该是可以通过某种手段去除的，但目前还没有查到，评论区知道的说一下哦。
下面上代码，使用滑动窗口进行切分，重复率设为0.1.

import os
import gdal
import numpy as np

# print(np.__path__)

#  读取tif数据集
def readTif(fileName):
    dataset = gdal.Open(fileName)
    if dataset == None:
        p
    return dataset


#  保存tif文件函数
def writeTiff(im_data, im_geotrans, im_proj, path):
    #指定数据类型
    if 'int8' in im_data.dtype.name: #输入的im_data是一个numpy。array类型，dtype的名字来对应gdal的datatype类型
        datatype = gdal.GDT_Byte
    elif 'int16' in im_data.dtype.name:
        datatype = gdal.GDT_UInt16
    else:
        datatype = gdal.GDT_Float32
    if len(im_data.shape) == 3: #形状 gadal读取出来的数据默认按波段、高、宽的形状排列[[每个波段对应的图像位置和值][][]]
        im_bands, im_height, im_width = im_data.shape
    elif len(im_data.shape) == 2: #如果只有二维意味着灰度，即[]，此时在外面再加一维表示波段即[[]]，类型为array，才可以识别为图片数组
        im_data = np.array([im_data]) #再加一维是可以直接实现的
        im_bands, im_height, im_width = im_data.shape
    # 创建文件
    driver = gdal.GetDriverByName("GTiff") #创建一个带有名字指代驱动器对象
    dataset = driver.Create(path, int(im_width), int(im_height), int(im_bands), datatype) #给出文件的路径，宽度，高度，波段数（倒过来）和数据类型，创建空文件，并确定开辟多大内存，像素值的类型用gdal类型
    #设置头文件信息：仿射变化参数，投影信息，数据体
    if (dataset != None):
        dataset.SetGeoTransform(im_geotrans)  # driver对象创建的dataset对象具有SetGeoTransform方法，写入仿射变换参数GEOTRANS
        dataset.SetProjection(im_proj)  # SetProjection写入投影im_proj
    for i in range(im_bands): #对每个波段进行处理 写入数据体
        dataset.GetRasterBand(i + 1).WriteArray(im_data[i]) #i从0开始，因此GetRasterBand(i+1),每个波段里用writearray写入图像信息数组im_data
    del dataset #写完之后释放内存空间


'''
滑动窗口裁剪函数
TifPath 影像路径
SavePath 裁剪后保存目录
CropSize 裁剪尺寸
RepetitionRate 重复率
'''


def TifCrop(TifPath, SavePath, CropSize, RepetitionRate): #tif文件裁剪函数，输入tif文件路径，保存路径，裁剪尺寸和重复率
    dataset_img = readTif(TifPath) #首先读取文件，采用gdal.open的方法
    width = dataset_img.RasterXSize #宽度是x方向上的size
    height = dataset_img.RasterYSize #高度是y方向上的size
    proj = dataset_img.GetProjection() #得到数据集的投影信息
    geotrans = dataset_img.GetGeoTransform() #得到数据集的地理仿射信息,是一个包含六个元素的元组

    img = dataset_img.ReadAsArray(0,0, width, height)  # 获取数据 Readasarray支持按块读取栅格，可以针对dataset也可以针对波段band,xoff = 0，yoff= 0xoff与yoff指定想要读取的部分原点位置在整张图像中距离全图的原点位置,这样就是读取全部数组信息了,GDAL读取的原点为左上角位置

    #  获取当前文件夹的文件个数len,并以len+1命名即将裁剪得到的图像
    new_name = len(os.listdir(SavePath)) + 1 #os.listdir(path)将path里面的文件夹列出来，os是operating system

    #  裁剪图片,重复率为RepetitionRate
    for i in range(int((height - CropSize * RepetitionRate) / (CropSize * (1 - RepetitionRate)))): #有重复率的裁减时，这个公式判断按照总高度height一共裁剪出多少幅图像 int函数将结果向下取整
        for j in range(int((width - CropSize * RepetitionRate) / (CropSize * (1 - RepetitionRate)))): #有重复率的裁剪时，按宽度一共裁剪出多少幅图像
            #  如果图像是单波段
            if (len(img.shape) == 2):
                cropped = img[
                          int(i * CropSize * (1 - RepetitionRate)): int(i * CropSize * (1 - RepetitionRate)) + CropSize,
                          int(j * CropSize * (1 - RepetitionRate)): int(j * CropSize * (1 - RepetitionRate)) + CropSize] #img[x:x+cropsize, y:y+cropsize] 多维数组的裁剪可以直接这么切片，
            # 如果图像是多波段
            else:
                cropped = img[:,
                          int(i * CropSize * (1 - RepetitionRate)): int(i * CropSize * (1 - RepetitionRate)) + CropSize,
                          int(j * CropSize * (1 - RepetitionRate)): int(j * CropSize * (1 - RepetitionRate)) + CropSize] #多波段，[波段，高，宽]的形状shape
            # 更新裁剪后的位置信息
            local_geotrans = list(geotrans) #每一张裁剪图的本地放射变化参数，0，3代表左上角坐标
            local_geotrans[0] = geotrans[0] + int(j * CropSize * (1 - RepetitionRate)) * geotrans[1]#分别更新为裁剪后的每一张局部图的左上角坐标，为滑动过的像素数量乘以分辨率
            local_geotrans[3] = geotrans[3] + int(i * CropSize * (1 - RepetitionRate)) * geotrans[5]
            local_geotrans = tuple(local_geotrans)
            # 写图像
            writeTiff(cropped, local_geotrans, proj, SavePath + "/%d.tif" % new_name) #数组、仿射变化参数、投影、保存路径
            #  文件名 + 1
            new_name = new_name + 1 #梅循环一次文件名的尾数+1

    # 向前裁剪最后一列
    for i in range(int((height - CropSize * RepetitionRate) / (CropSize * (1 - RepetitionRate)))):
        if (len(img.shape) == 2):
            cropped = img[int(i * CropSize * (1 - RepetitionRate)): int(i * CropSize * (1 - RepetitionRate)) + CropSize,
                      (width - CropSize): width] #刚好一个cropsize的宽，也就是最后一列跟前面重复率可能为repetition~1
        else:
            cropped = img[:,
                      int(i * CropSize * (1 - RepetitionRate)): int(i * CropSize * (1 - RepetitionRate)) + CropSize,
                      (width - CropSize): width]
        # 更新裁剪后的位置信息
        local_geotrans = list(geotrans)  # 每一张裁剪图的本地放射变化参数，0，3代表左上角坐标
        local_geotrans[0] = geotrans[0] + (width - CropSize) * geotrans[1]  # 分别更新为裁剪后的每一张局部图的左上角坐标
        local_geotrans[3] = geotrans[3] + int(i * CropSize * (1 - RepetitionRate)) * geotrans[5]
        local_geotrans = tuple(local_geotrans)
        # 写图像
        writeTiff(cropped, local_geotrans, proj, SavePath + "/%d.tif" % new_name)
        new_name = new_name + 1

    # 向前裁剪最后一行
    for j in range(int((width - CropSize * RepetitionRate) / (CropSize * (1 - RepetitionRate)))):
        if (len(img.shape) == 2):
            cropped = img[(height - CropSize): height,
                      int(j * CropSize * (1 - RepetitionRate)): int(j * CropSize * (1 - RepetitionRate)) + CropSize]#刚好一个cropsize的高，也就是最后一列跟前面重复率可能为repetition~1
        else:
            cropped = img[:,
                      (height - CropSize): height,
                      int(j * CropSize * (1 - RepetitionRate)): int(j * CropSize * (1 - RepetitionRate)) + CropSize]
        # 更新裁剪后的位置信息
        local_geotrans = list(geotrans)  # 每一张裁剪图的本地放射变化参数，0，3代表左上角坐标
        local_geotrans[0] = geotrans[0] + int(j * CropSize * (1 - RepetitionRate)) * geotrans[1]  # 分别更新为裁剪后的每一张局部图的左上角坐标
        local_geotrans[3] = geotrans[3] + (height - CropSize) * geotrans[5]
        local_geotrans = tuple(local_geotrans)
        #写文件
        writeTiff(cropped, local_geotrans, proj, SavePath + "/%d.tif" % new_name)
        #  文件名 + 1
        new_name = new_name + 1

    # 裁剪右下角
    if (len(img.shape) == 2):
        cropped = img[(height - CropSize): height,
                  (width - CropSize): width] #保证右下角也完全裁剪到，因为这种滑动窗口类似有向裁剪，在终点处的部分不一定能完全被裁剪到，所以右下角是各个位置为起点的终点，需要单独考虑
    else:
        cropped = img[:,
                  (height - CropSize): height,
                  (width - CropSize): width]
    # 更新裁剪后的位置信息
    local_geotrans = list(geotrans)  # 每一张裁剪图的本地放射变化参数，0，3代表左上角坐标
    local_geotrans[0] = geotrans[0] +  (width - CropSize) * geotrans[1]  # 分别更新为裁剪后的每一张局部图的左上角坐标
    local_geotrans[3] = geotrans[3] + (height - CropSize) * geotrans[5]
    local_geotrans = tuple(local_geotrans)
    writeTiff(cropped, local_geotrans, proj, SavePath + "/%d.tif" % new_name)
    new_name = new_name + 1
#  将影像1裁剪为重复率为0.1的256×256的数据集
TifCrop(r"D:\dataset\五布河\DOMDSM\五布河.TIF", r"D:\deep_learning\Dataset\RSdata\Wubu", 2560, 0.1)