python 遥感图像融合 python遥感影像分类

主要功能

批量读取文件，借助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错误

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笔记】分块下载——像元行列数