【Python笔记】遥感影像(tif)批量分块

• 主要功能

批量读取文件，借助GDAL以及numpy分块遥感数据——固定行列像元数量

• 环境配置

主要版本如下：

python                   3.9.7

gdal                       3.4.1

numpy                   1.21.5

pyproj                    3.3.0

• 代码实现

1、数据载入

a、载入相关包

        os.environ添加环境变量是因为碰到错误，正常可不需要GDAL找不到proj.db错误https://blog.csdn.net/weixin_40694662/article/details/124183483

b、设定影像所在文件夹的绝对路径 tifDir，结果输出路径 outPath

 

 

c、读取文件夹内后缀为 .tif的文件，构成文件名数组tifs

d、for循环去掉 tifs 元素后四位（.tif）,构成日期列表datelist

from osgeo import gdal
import numpy as np
import os
#os.environ['PROJ_LIB'] = r'E:\ProgramData\Anaconda3\envs\qt39\Lib\site-packages\pyproj\proj_dir\share\proj'
# os.environ['GDAL_DATA'] = r'E:\ProgramData\Anaconda3\envs\qt39\Library\share'
import random
from tqdm import tqdm


# 分块影像所在文件夹，不能有中文
tifDir = r"E:\pyimg\tif2csv\S2SR10mallband3"
# 输出的文件夹，不能有中文，如果文件夹不存在则会被创建
outPath = r"E:\pyimg\tif2csv\S2SR10mallband3tile"
if not os.path.exists(outPath):
    os.makedirs(outPath)

tifs = [i for i in os.listdir(tifDir) if i.endswith(".tif")]

print("有 %s 个tif文件" % len(tifs))
print("tifs",tifs)
datelist1 = []
for i in tifs:
    datelist1.append(i[:-4])
datelist = list(set(datelist1))
datelist.sort(key=datelist1.index)
print("有 %s 个日期" % len(datelist))
print("datelist" , datelist)

 2、数据分块

a、设置分块大小，这里分块为正方形，行列像元数目size设置为256

b、for循环遍历文件，gdal打开，获取相关元数据，以及数据矩阵im_data(波段，行数，列数)

c、嵌套for循环可分割的行数、列数，1）分割数据矩阵im_data 2）设置输出文件名为“行数-列数”3）计算新的原点坐标，复制原图的投影，连同分割的im_data写入新文件

# 定义切图的大小（矩形框）
size = 256

for img in range(len(datelist)):
    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
    im_data = in_ds.ReadAsArray()  # 获取数据


    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 = i * size
            offset_y = j * size
            ## 从每个波段中切需要的矩形框内的数据(注意读取的矩形框不能超过原图大小)
            b_ysize = min(width - offset_y, size)
            b_xsize = min(height - offset_x, size)

            # print("width:%d     height:%d    offset_x:%d    offset_y:%d     b_xsize:%d     b_ysize:%d" % (width, height, offset_x, offset_y, b_xsize, b_ysize))
            # print(im_data.shape)
            #裁剪影像矩阵
            out_allband = im_data[:,offset_x:offset_x+b_xsize,offset_y:offset_y+b_ysize]
            # print(out_allband.shape)

            # 获取Tif的驱动，为创建切出来的图文件做准备
            gtif_driver = gdal.GetDriverByName("GTiff")
            file = outPath+"\\"+tifs[img][:-4]+"-"+str(offset_x).zfill(10)+"-"+str(offset_y).zfill(10)+".tif"

            # 创建切出来的要存的文件
            out_ds = gtif_driver.Create(file, b_ysize, b_xsize, outbandsize, datatype)
            # print("create new tif file succeed")

            # 获取原图的原点坐标信息
            ori_transform = in_ds.GetGeoTransform()
            # if ori_transform:
                # print(ori_transform)
                # print("Origin = ({}, {})".format(ori_transform[0], ori_transform[3]))
                # print("Pixel Size = ({}, {})".format(ori_transform[1], ori_transform[5]))

            # 读取原图仿射变换参数值
            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_y * w_e_pixel_resolution
            top_left_y = top_left_y + offset_x * 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):
                out_ds.GetRasterBand(ii + 1).WriteArray(out_allband[ii])

            # 将缓存写入磁盘
            out_ds.FlushCache()
            # print("FlushCache succeed")
            del out_ds

d、结果 输出的数据结果命名方式与GEE分块下载相同,GEE数据详情可参阅【GEE笔记】分块下载——像元行列数

​