随着科技的进步，遥感影像包含的信息越来越多，存储空间也变得很大，这就导致我们在处理影像时软件会非常的卡。同时目前很火的深度学习也需要对影像分割后制作样本集，所以今天给大家分享下如何使用Python的GDAL库对遥感影像进行分幅裁剪！
一、导入需要的三方库
tkinter.filedialog是为了使用窗口打开影像，这样就不用每次使用时都修改一下影像路径了。numpy库是为了以数组的形式读取波段。
import os
from osgeo import gdal
import numpy as np
import tkinter.filedialog
二、读取影像的基本信息
def Get_data(filepath):
"""
:param filepath: 输入数据路径
:return: 输出影像的基本信息
"""
ds = gdal.Open(filepath) # 打开数据集dataset
ds_width = ds.RasterXSize # 获取数据宽度
ds_height = ds.RasterYSize # 获取数据高度
ds_bands = ds.RasterCount # 获取波段数
ds_geo = ds.GetGeoTransform() # 获取仿射地理变换参数
ds_prj = ds.GetProjection() # 获取投影信息
print("影像的宽度为：" + str(ds_width))
print("影像的高度为：" + str(ds_height))
print("影像的波段数为：" + str(ds_bands))
print("仿射地理变换参数为：" + str(ds_geo))
print("投影坐标系为：" + str(ds_prj))
# data = ds.ReadAsArray(0, 0, ds_width, ds_height) # 以数组的形式读取整个数据集三、分幅裁剪
1.计算裁剪的宽度、高度
raw是多少行，col是多少列。
raw_frame = int(ds_height / raw)
# 计算每幅图像的高度
col_frame = int(ds_width / col)
# 计算每幅图像的宽度2.计算分幅左上角的像素坐标
j,k是通过for循环每一行每一列，这个在后面的完整代码中会有展示
left_x = j * raw_frame
left_y = k * col_frame
# 计算当前影像的左上角像素坐标3.获取新的投影和仿射地理变换参数
ds_geo = ds.GetGeoTransform() # 获取仿射地理变换参数
top_left_x = ds_geo[0] # 原始影像左上角x的投影坐标
top_left_y = ds_geo[3] # 原始影像左上角y投影坐标
top_left_x = top_left_x + left_y * ds_geo[1]
top_left_y = top_left_y + left_x * ds_geo[5]
# 计算得到当前影像的左上角投影坐标
ds_geo = (top_left_x, ds_geo[1], ds_geo[2], top_left_y, ds_geo[4], ds_geo[5])
# 新影像的仿射地理变换参数
ds_result.SetGeoTransform(ds_geo) # 导入仿射地理变换参数
ds_result.SetProjection(ds.GetProjection()) # 导入投影信息4.创建新的tif循环写入各个波段
driver = gdal.GetDriverByName('GTiff') # 载入数据驱动，用于存储内存中的数组
ds_result = driver.Create(out_path+"%s_%s.tif" % (j+1, k+1), col_frame, raw_frame,
bands=ds_bands, eType=gdal.GDT_Float64) # 创建空tif
for i in range(1, ds_bands+1):
array_band = ds.GetRasterBand(i).ReadAsArray(left_y, left_x, col_frame,raw_frame).astype(np.float64)
# 根据左上角的像素坐标和幅宽读取指定区域内的数据
ds_result.GetRasterBand(i).SetNoDataValue(0) # 将无效值设为0
ds_result.GetRasterBand(i).WriteArray(array_band) # 将每个波段写入新的文件中5.删除空值
影像一般都不是正正好好的矩形，所以裁剪时会有很多没有数据的影像，所以我们要把这些影像删除。
if array_band.any() == 0:
ds_result = None
print("此幅影像为空值，已删除！")
os.remove(out_path+"%s_%s.tif" % (j+1, k+1))四、完整代码
# -*- coding: utf-8 -*-
"""
@Time ： 2023/8/14 11:34
@Auth ： RS迷途小书童
@File ：Image Framing and Cropping.py
@IDE ：PyCharm
@Purpose ：遥感影像分幅裁剪
"""
import os
from osgeo import gdal
import numpy as np
import tkinter.filedialog
def Get_data(filepath):
"""
:param filepath: 输入数据路径
:return: 输出影像的基本信息
"""
ds = gdal.Open(filepath) # 打开数据集dataset
ds_width = ds.RasterXSize # 获取数据宽度
ds_height = ds.RasterYSize # 获取数据高度
ds_bands = ds.RasterCount # 获取波段数
ds_geo = ds.GetGeoTransform() # 获取仿射地理变换参数
ds_prj = ds.GetProjection() # 获取投影信息
print("影像的宽度为：" + str(ds_width))
print("影像的高度为：" + str(ds_height))
print("影像的波段数为：" + str(ds_bands))
print("仿射地理变换参数为：" + str(ds_geo))
print("投影坐标系为：" + str(ds_prj))
# data = ds.ReadAsArray(0, 0, ds_width, ds_height) # 以数组的形式读取整个数据集
def Frame_image(filepath, out_path, raw, col):
"""
:param filepath: 输入需要分幅裁剪的影像路径
:param out_path: 输出文件夹即可，名称固定为行列
:param raw: 需要分幅的行数
:param col: 需要分幅的列数
:return: None
"""
ds = gdal.Open(filepath) # 打开数据集dataset
ds_width = ds.RasterXSize # 获取数据宽度
ds_height = ds.RasterYSize # 获取数据高度
ds_bands = ds.RasterCount # 获取波段数
for j in range(0, raw):
for k in range(0, col):
print("正在裁剪第%s行第%s列......" % (j+1, k+1))
raw_frame = int(ds_height / raw)
# 计算每幅图像的高度
col_frame = int(ds_width / col)
# 计算每幅图像的宽度
left_x = j * raw_frame
left_y = k * col_frame
# 计算当前影像的左上角像素坐标
raw_frame = min(ds_height-left_x, raw_frame)
col_frame = min(ds_width-left_y, col_frame)
# 防止幅宽超过整幅影像
driver = gdal.GetDriverByName('GTiff') # 载入数据驱动，用于存储内存中的数组
ds_result = driver.Create(out_path+"%s_%s.tif" % (j+1, k+1), col_frame, raw_frame,
bands=ds_bands, eType=gdal.GDT_Float64) # 创建空tif
ds_geo = ds.GetGeoTransform() # 获取仿射地理变换参数
top_left_x = ds_geo[0] # 原始影像左上角x的投影坐标
top_left_y = ds_geo[3] # 原始影像左上角y投影坐标
top_left_x = top_left_x + left_y * ds_geo[1]
top_left_y = top_left_y + left_x * ds_geo[5]
# 计算得到当前影像的左上角投影坐标
ds_geo = (top_left_x, ds_geo[1], ds_geo[2], top_left_y, ds_geo[4], ds_geo[5])
# 新影像的仿射地理变换参数
ds_result.SetGeoTransform(ds_geo) # 导入仿射地理变换参数
ds_result.SetProjection(ds.GetProjection()) # 导入投影信息
array_band = []
for i in range(1, ds_bands+1):
array_band = ds.GetRasterBand(i).ReadAsArray(left_y, left_x, col_frame, raw_frame).astype(np.float64)
# 根据左上角的像素坐标和幅宽读取指定区域内的数据
ds_result.GetRasterBand(i).SetNoDataValue(0) # 将无效值设为0
ds_result.GetRasterBand(i).WriteArray(array_band) # 将每个波段写入新的文件中
if array_band.any() == 0:
ds_result = None
print("此幅影像为空值，已删除！")
os.remove(out_path+"%s_%s.tif" % (j+1, k+1))
ds_result = None
if __name__ == "__main__":
fn_input = open(tkinter.filedialog.askopenfilename(title='选择图片', filetypes=[('所有文件', '.*'), ('JPG', '.jpg'), ('JPG', '.jpeg'), ('TIFF', '.tif'), ('DAT', '.dat'), ('PNG', '.png')]), 'rb')
path = 'B:/Personal/Vegetation enhancement/layer'
out_path1 = 'B:/Personal/Python_分幅/'
# 输出文件夹，不需要文件名，通过行列号命名
raw1 = int(input("请输入行数："))
col1 = int(input("请输入列数："))
Frame_image(fn_input.name, out_path1, raw1, col1)
今天主要分享的是遥感影像的分幅裁剪，大家可以用这段代码减少数据量，也可以用它制作样本集。如果大家在学习Python或者RS时有什么问题，可以随时留言交流！如果大家对批量处理有兴趣同样可以留言给博主，博主会分享相关代码以供学习！