上一篇文章简单介绍了波段的叠加，本文对叠加的波段进行影像裁剪，输出是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')
本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。
