用gdal分块读取和裁剪卫星影像

起因

最近需要将高分-2的高分辨率影像裁剪成小块用于制作深度学习的训练样本，因而研究了一下如何使用python 代码来批量对影像进行切块。

这里主要需要使用到的是gdal和numpy两个包

通过如下代码可以安装这两个包：

pip install gdal
pip install numpy
#or
conda install gdal
conda install numpy

导入相应的库

from osgeo import gdal
import numpy as np
import os
import random
from tqdm import tqdm

获取待处理的影像的路径

#step1: read the image and set the outpath

# 分块影像所在文件夹，不能有中文
tifDir = r"E:\GF-image"
# 输出的文件夹，不能有中文，如果文件夹不存在则会被创建
outPath = r"E:\GF-imageout"
if not os.path.exists(outPath):
    os.makedirs(outPath)
#将tif影像都找出来得到一个列表 
tifs = [i for i in os.listdir(tifDir) if i.endswith(".tif")]
 
print("有 %s 个tif文件" % len(tifs))
print("tifs",tifs)

边读取边导出

我这里设定导出影像大小为256*256,故size = 256.

# step2: clip the image to 256x256
# 定义切图的大小（矩形框）
size = 256
for img in range(len(tifs)):
    # print("正在分割：",tifs[img])
    in_ds = gdal.Open(tifDir + "\\" + tifs[img])              # 读取要切的原图
 
    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
    dataTYpe = gdal.GetDataTypeName(datatype)
 
    col_num = int(width / size)  # 宽度可以分成几块
    row_num = int(height / size)  # 高度可以分成几块
    if (width % size != 0):
        col_num += 1
    if (height % size != 0):
        row_num += 1
 
    print("row_num:%d   col_num:%d" % (row_num, col_num))
    for i in tqdm(range(row_num)):  
        for j in range(col_num):
            offset_x = j * size
            offset_y = i * size
            ## 从每个波段中切需要的矩形框内的数据(注意读取的矩形框不能超过原图大小)
            b_xsize = min(width - offset_x, size)
            b_ysize = min(height - offset_y, size)
 
            #裁剪影像矩阵
            im_data = in_ds.ReadAsArray(offset_x, offset_y, b_xsize, b_ysize)
            out_allband = im_data[:3,:,:]
 
            # 获取Tif的驱动，为创建切出来的图文件做准备
            gtif_driver = gdal.GetDriverByName("GTiff")
            file = outPath+"\\"+str(i)+"-"+str(j)+".tif"
 
            # 创建切出来的要存的文件
            out_ds = gtif_driver.Create(file, b_xsize, b_ysize, 3, gdal.GDT_UInt16)
            # print("create new tif file succeed")
 
            # 获取原图的原点坐标信息
            ori_transform = in_ds.GetGeoTransform()
 
            # 读取原图仿射变换参数值
            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 + offset_x * w_e_pixel_resolution
            top_left_y = top_left_y + offset_y * 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())
 
            # 写入目标文件
            for ii in range(outbandsize-1):
                out_ds.GetRasterBand(ii+1).WriteArray(out_allband[ii])
 
            # 将缓存写入磁盘
            out_ds.FlushCache()
            # print("FlushCache succeed")
            del out_ds

1. 这里面需要注意的是gdal是以width轴为x轴，以height轴为y轴。读取完一个block以后，对应的左上角点坐标（xoffset,yoffset）需要更新，（xsize,ysize）要和它的x,y保持对应，同时需要更新导出块的仿射变换信息。
2. gdal读取波段是从1开始的
3. 每个block导出以后需要清空缓存，否则内存将爆炸。