【Python&RS】基于GDAL遥感影像分幅裁剪

标签：RS Python frame col GDAL geo ds 影像 left

随着科技的进步，遥感影像包含的信息越来越多，存储空间也变得很大，这就导致我们在处理影像时软件会非常的卡。同时目前很火的深度学习也需要对影像分割后制作样本集，所以今天给大家分享下如何使用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时有什么问题，可以随时留言交流！如果大家对批量处理有兴趣同样可以留言给博主，博主会分享相关代码以供学习！

标签：RS,Python,frame,col,GDAL,geo,ds,影像,left
From： https://www.cnblogs.com/RSran/p/17631881.html