python影像裁剪并保存成tiff格式（规则网格法）

上一篇文章简单介绍了波段的叠加，本文对叠加的波段进行影像裁剪，输出是channels×256×256的tiff格式。本文构建了相应的mask矩阵加在了channel的最后一个通道。

数据准备

波段的读取可以参考python批量读取landsat8的波段，具体的函数整理数据的函数在这篇文章都介绍的很详细，这里不再重复，我们直接用它的返回列表。

波段叠加的具体方式参考python矩阵堆叠-实现遥感影像波段叠加，这些都是博主刚刚更新过的内容，这里直接拿来用。

from read_landsat8 import read_landsat8_bands
import numpy as np
from osgeo import gdal_array

def get_img(base_path):
    """
    叠加波段的简单demo
    与mask矩阵一起构建成3维数组结构"""
    bands = read_landsat8_bands(base_path)

    # 读取数据
    B1_gdal = gdal_array.LoadFile(bands[0][0])
    B2_gdal = gdal_array.LoadFile(bands[0][1])
    B3_gdal = gdal_array.LoadFile(bands[0][2])

    # 转化成ndarray形式
    B1_np = np.array(B1_gdal)
    B2_np = np.array(B2_gdal)
    B3_np = np.array(B3_gdal)
    print(B1_np.shape)

    B123 = np.stack([B1_np, B2_np, B3_np], axis=0)
    print(B123.shape)  # 3,7301,7341

    # 构建0-1 mask矩阵
    height = B123.shape[1]
    width = B123.shape[2]
    mask = np.random.randint(0, 2, (1, height, width))

    # 按照通道堆叠
    img = np.concatenate([B123, mask], axis=0)

    return img

这里得到的img其实是一个3维数组，在本例中，它的维度是
4×7301×7341

envi查看一下堆叠后的图像

图像裁剪

有了上面构建的img数据，再设定一个裁剪尺寸，我们就可以进行图像裁剪了。

def crop_img(img, cropsize):
    """
    裁剪图像为指定格式并保存成tiff
    输入为array形式的数组
    """
    num = 0
    height = img.shape[1]
    width = img.shape[2]

    # 从左上开始裁剪
    for i in range(int((height) / (cropsize))):  # 行裁剪次数
        for j in range(int((width) / (cropsize))):  # 列裁剪次数
            cropped = img[:,  # 通道不裁剪
                      i * cropsize: i * cropsize + cropsize,
                      j * cropsize: j * cropsize + cropsize,
                      ]

            num = num + 1
            target = 'tiff_crop' + '/cropped{n}.tif'.format(n=num)
            out = gdal_array.SaveArray(cropped, target, format="GTiff")

    # #  向前裁剪最后一列
    for i in range(int((height) / (cropsize))):
        for j in range(int((width) / (cropsize))):
            cropped = img[:,  # 通道不裁剪
                      i * cropsize: i * cropsize + cropsize,  # 所有行
                      width - cropsize: width,  # 最后256列
                      ]

            num = num + 1
            target = 'tiff_crop' + '/cropped{n}.tif'.format(n=num)
            out = gdal_array.SaveArray(cropped, target, format="GTiff")

    # #  向前裁剪最后一行
    for i in range(int((height) / (cropsize))):
        for j in range(int((width) / (cropsize))):
            cropped = img[:,  # 通道不裁剪
                      height - cropsize: height,  # 最后256行
                      j * cropsize: j * cropsize + cropsize,  # 所有列
                      ]

            num = num + 1
            target = 'tiff_crop' + '/cropped{n}.tif'.format(n=num)
            out = gdal_array.SaveArray(cropped, target, format="GTiff")
   
    # 裁剪右下角
    cropped = img[:,  # 通道不裁剪
              height - cropsize: height,
              width - cropsize: width,
              ]

    num = num + 1
    target = 'tiff_crop' + '/cropped{n}.tif'.format(n=num)
    gdal_array.SaveArray(cropped, target, format="GTiff")

这里img就是前面得到的3维数组，cropsize是裁剪尺寸。

在裁剪时没有丢弃内容，而是将所有行列都裁剪到了。从左上角开始裁剪，最后右下角的边角料也有专门的处理。

裁剪后保存成了tiff格式。

测试

设置裁剪尺寸为256，看一下输出结果

if __name__ == '__main__':
    base_path = 'data'
    img = get_img(base_path)

    cropsize = 256
    crop_img(img, cropsize)
    print('finish')

out：

(7301, 7341)
(3, 7301, 7341)
finish

我们来到保存裁剪后图像的目标文件夹：
这张tif比较大，一共裁出了2000多个小tif

用envi打开其中一张康康

所有代码

from read_landsat8 import read_landsat8_bands
import numpy as np
from osgeo import gdal_array


def crop_img(img, cropsize):
    """
    裁剪图像为指定格式并保存成tiff
    输入为array形式的数组
    """
    num = 0
    height = img.shape[1]
    width = img.shape[2]
    print(height)
    print(width)

    # 从左上开始裁剪
    for i in range(int((height) / (cropsize))):  # 行裁剪次数
        for j in range(int((width) / (cropsize))):  # 列裁剪次数
            cropped = img[:,  # 通道不裁剪
                      i * cropsize: i * cropsize + cropsize,
                      j * cropsize: j * cropsize + cropsize,
                      ]

            num = num + 1
            target = 'tiff_crop' + '/cropped{n}.tif'.format(n=num)
            gdal_array.SaveArray(cropped, target, format="GTiff")

    # #  向前裁剪最后一列
    for i in range(int((height) / (cropsize))):
        for j in range(int((width) / (cropsize))):
            cropped = img[:,  # 通道不裁剪
                      i * cropsize: i * cropsize + cropsize,  # 所有行
                      width - cropsize: width,  # 最后256列
                      ]

            num = num + 1
            target = 'tiff_crop' + '/cropped{n}.tif'.format(n=num)
            gdal_array.SaveArray(cropped, target, format="GTiff")

    # #  向前裁剪最后一行
    for i in range(int((height) / (cropsize))):
        for j in range(int((width) / (cropsize))):
            cropped = img[:,  # 通道不裁剪
                      height - cropsize: height,  # 最后256行
                      j * cropsize: j * cropsize + cropsize,  # 所有列
                      ]

            num = num + 1
            target = 'tiff_crop' + '/cropped{n}.tif'.format(n=num)
            gdal_array.SaveArray(cropped, target, format="GTiff")

    # 裁剪右下角
    cropped = img[:,  # 通道不裁剪
              height - cropsize: height,
              width - cropsize: width,
              ]

    num = num + 1
    target = 'tiff_crop' + '/cropped{n}.tif'.format(n=num)
    gdal_array.SaveArray(cropped, target, format="GTiff")



def get_img(base_path):
    """
    叠加波段的简单demo
    与mask矩阵一起构建成3维数组结构"""
    bands = read_landsat8_bands(base_path)

    # 读取数据
    B1_gdal = gdal_array.LoadFile(bands[0][0])
    B2_gdal = gdal_array.LoadFile(bands[0][1])
    B3_gdal = gdal_array.LoadFile(bands[0][2])

    # 转化成ndarray形式
    B1_np = np.array(B1_gdal)
    B2_np = np.array(B2_gdal)
    B3_np = np.array(B3_gdal)
    print(B1_np.shape)

    B123 = np.stack([B1_np, B2_np, B3_np], axis=0)
    print(B123.shape)  # 3,7301,7341

    # 构建0-1 mask矩阵
    height = B123.shape[1]
    width = B123.shape[2]
    mask = np.random.randint(0, 2, (1, height, width))

    # 按照通道堆叠
    img = np.concatenate([B123, mask], axis=0)

    return img



if __name__ == '__main__':
    base_path = 'data'

    img = get_img(base_path)

    cropsize = 256
    crop_img(img, cropsize)
    print('finish')