用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)
  • 14
    点赞
  • 37
    收藏
    觉得还不错? 一键收藏
  • 10
    评论
GDAL(Geospatial Data Abstraction Library)是一个开源的地理空间数据读写库。它可以读取和处理多种格式的地理空间数据,包括TIFF图像。TIFF(Tagged Image File Format)是一种常用的图像格式,它支持多种数据类型、多通道和元数据。GDAL通过使用TIFF库来读取TIFF文件中的数据。 在使用GDAL读取TIFF文件时,我们需要先安装GDAL库及相关依赖库。安装完成后,可以使用代码来读取TIFF文件中的数据。 使用GDAL库读取TIFF文件的基本步骤如下: 1.导入GDAL库和相关依赖库 import gdal 2.打开TIFF文件 ds = gdal.Open('filename.tif') 其中,'filename.tif'是需要读取的TIFF文件名。 3.获取TIFF文件的元数据 获取TIFF文件的元数据,包括文件的坐标系、分辨率、波段数等。 projection = ds.GetProjection() # 坐标系信息 geo_transform = ds.GetGeoTransform() # 分辨率等信息 band_nums = ds.RasterCount # 波段数 4.读取TIFF文件中的数据 读取TIFF文件中的数据,可以使用多种方式,包括读取整个文件,读取指定区域、指定波段等。例如,读取第一波段的整个数据: data = ds.GetRasterBand(1).ReadAsArray() 其中,data是一个二维数组,包含了TIFF文件中第一波段的全部数据。 5.关闭TIFF文件 使用完TIFF文件后记得关闭。 ds.Close() 以上就是使用GDAL库读取TIFF文件的基本步骤。在实际应用中,还可以对TIFF文件进行裁剪、投影变换等操作,以满足不同的需求。
评论 10
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值