1、GDAL进行图像一块区域的裁剪和大图像裁剪为多个小图像

1 代码实现(一小块儿区域裁剪)

# 将一幅遥感影像拆建为多张
import gdal

# 读取要切的原图
in_ds = gdal.Open("zx1.jpg")
print("open tif file succeed")
width = in_ds.RasterXSize  # 获取数据宽度
height = in_ds.RasterYSize  # 获取数据高度
outbandsize = in_ds.RasterCount  # 获取数据波段数
im_geotrans = in_ds.GetGeoTransform()  # 获取仿射矩阵信息
im_proj = in_ds.GetProjection()  # 获取投影信息
datatype = in_ds.GetRasterBand(1).DataType
im_data = in_ds.ReadAsArray()  # 获取数据

# 读取原图中的每个波段
in_band1 = in_ds.GetRasterBand(1)
in_band2 = in_ds.GetRasterBand(2)
in_band3 = in_ds.GetRasterBand(3)

# 定义切图的起始点坐标
offset_x = 0
offset_y = 0

# offset_x = width/2  # 这里是随便选取的,可根据自己的实际需要设置
# offset_y = height/2

# 定义切图的大小(矩形框)
block_xsize = 60  # 行
block_ysize = 60  # 列

k = 0
for i in range(width // block_xsize):
    for j in range(height // block_xsize):
        out_band1 = in_band1.ReadAsArray(i * block_xsize, j * block_xsize, block_xsize, block_ysize)
        out_band2 = in_band2.ReadAsArray(i * block_xsize, j * block_xsize, block_xsize, block_ysize)
        out_band3 = in_band3.ReadAsArray(i * block_xsize, j * block_xsize, block_xsize, block_ysize)
        print(out_band3)
        k += 1

        ## 从每个波段中切需要的矩形框内的数据(注意读取的矩形框不能超过原图大小)
        # out_band1 = in_band1.ReadAsArray(offset_x, offset_y, block_xsize, block_ysize)
        # out_band2 = in_band2.ReadAsArray(offset_x, offset_y, block_xsize, block_ysize)
        # out_band3 = in_band3.ReadAsArray(offset_x, offset_y, block_xsize, block_ysize)

        # 获取Tif的驱动,为创建切出来的图文件做准备
        gtif_driver = gdal.GetDriverByName("GTiff")

        # 创建切出来的要存的文件(3代表3个不都按,最后一个参数为数据类型,跟原文件一致)
        out_ds = gtif_driver.Create(str(k) + 'clip4.tif', block_xsize, block_ysize, outbandsize, datatype)
        # print("create new tif file succeed")

        # 获取原图的原点坐标信息,# 获取仿射矩阵信息
        ori_transform = in_ds.GetGeoTransform()
        if ori_transform:
            print(ori_transform)
            print("Origin = ({}, {})".format(ori_transform[0], ori_transform[3]))
            print("Pixel Size = ({}, {})".format(ori_transform[1], ori_transform[5]))

        # 读取原图仿射变换参数值
        top_left_x = ori_transform[0]  # 左上角x坐标
        w_e_pixel_resolution = ori_transform[1]  # 东西方向像素分辨率
        top_left_y = ori_transform[3]  # 左上角y坐标
        n_s_pixel_resolution = ori_transform[5]  # 南北方向像素分辨率

        # 根据反射变换参数计算新图的原点坐标
        top_left_x = top_left_x + i * block_xsize * w_e_pixel_resolution
        top_left_y = top_left_y + j * block_xsize * n_s_pixel_resolution

        # 将计算后的值组装为一个元组,以方便设置
        dst_transform = (top_left_x, ori_transform[1], ori_transform[2], top_left_y, ori_transform[4], ori_transform[5])

        # 设置裁剪出来图的原点坐标
        out_ds.SetGeoTransform(dst_transform)

        # 设置SRS属性(投影信息)
        out_ds.SetProjection(in_ds.GetProjection())

        # 写入目标文件
        out_ds.GetRasterBand(1).WriteArray(out_band1)
        out_ds.GetRasterBand(2).WriteArray(out_band2)
        out_ds.GetRasterBand(3).WriteArray(out_band3)

        # 将缓存写入磁盘
        out_ds.FlushCache()
        print("FlushCache succeed")

        # 计算统计值
        # for i in range(1, 3):
        #     out_ds.GetRasterBand(i).ComputeStatistics(False)
        # print("ComputeStatistics succeed")

        del out_ds

        print("End!")

在这里插入图片描述

2 将一大块区域裁剪为小块图像(损失了坐标)

import os
import gdal
import numpy as np


#  读取tif数据集
def readTif(fileName):
    dataset = gdal.Open(fileName)
    if dataset == None:
        print(fileName + "文件无法打开")
    return dataset


#  保存tif文件函数
def writeTiff(im_data, im_geotrans, im_proj, path):
    if 'int8' in im_data.dtype.name:
        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:
        im_bands, im_height, im_width = im_data.shape
    elif len(im_data.shape) == 2:
        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)
    if (dataset != None):
        dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
        dataset.SetProjection(im_proj)  # 写入投影
    for i in range(im_bands):
        dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
    del dataset

def TifCrop(TifPath, SavePath, CropSize, RepetitionRate):
    '''
    滑动窗口裁剪函数
    TifPath 影像路径
    SavePath 裁剪后保存目录
    CropSize 裁剪尺寸
    RepetitionRate 重复率
    '''
    dataset_img = readTif(TifPath)
    width = dataset_img.RasterXSize
    height = dataset_img.RasterYSize
    proj = dataset_img.GetProjection()
    geotrans = dataset_img.GetGeoTransform()
    img = dataset_img.ReadAsArray(0, 0, width, height)  # 获取数据

    #  获取当前文件夹的文件个数len,并以len+1命名即将裁剪得到的图像
    #new_name = len(os.listdir(SavePath)) + 1
    new_name = 1
    #  裁剪图片,重复率为RepetitionRate

    for i in range(int((height - CropSize * RepetitionRate) / (CropSize * (1 - RepetitionRate)))):
        print("row_number:%d"%i)
        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]
            #  如果图像是多波段
            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]
            #  写图像
            writeTiff(cropped, geotrans, proj, SavePath + "/n%d.tif" % new_name)
            #  文件名 + 1
            new_name = new_name + 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]
        else:
            cropped = img[:,
                      int(i * CropSize * (1 - RepetitionRate)): int(i * CropSize * (1 - RepetitionRate)) + CropSize,
                      (width - CropSize): width]
        #  写图像
        writeTiff(cropped, geotrans, proj, SavePath + "/n%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]
        else:
            cropped = img[:,
                      (height - CropSize): height,
                      int(j * CropSize * (1 - RepetitionRate)): int(j * CropSize * (1 - RepetitionRate)) + CropSize]
        writeTiff(cropped, geotrans, proj, SavePath + "/n%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]
    writeTiff(cropped, geotrans, proj, SavePath + "/n%d.tif" % new_name)
    new_name = new_name + 1

if __name__ == '__main__':
    #  将影像1裁剪为重复率为0.1的256×256的数据集
    TifCrop(r"H:\01HTutorWork\3GF2\2DataAndLabel\4Gaofen2\1fujian_shuncang\2imageLabel\zhencaiseGF2_PMS2_E118.1_N27.0_20190919_L1A0004254305-PAN2_ORTHO_PSH1.tif",
            r"H:\01HTutorWork\3GF2\2DataAndLabel\4Gaofen2\1fujian_shuncang\5dataset16\1images", 16, 0)
    TifCrop(r"H:\01HTutorWork\3GF2\2DataAndLabel\4Gaofen2\1fujian_shuncang\2imageLabel\sc254305.tif",
            r"H:\01HTutorWork\3GF2\2DataAndLabel\4Gaofen2\1fujian_shuncang\5dataset16\2labels", 16, 0)

参考资料

遥感图像裁剪
https://blog.csdn.net/zsc201825/article/details/89359995?utm_medium=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-2.control&dist_request_id=1328641.12274.16155548391037017&depth_1-utm_source=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-2.control

gdal图像格式转换
https://blog.csdn.net/godenlove007/article/details/8864763

  • 5
    点赞
  • 26
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

晓码bigdata

如果文章给您带来帮助,感谢打赏

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

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

余额充值