一、问题引入

        最近在做遥感语义分割的实验时，发现样本的制作是一个比较棘手的问题，因此想要整理下代码，以供学习交流。

         我的原始遥感图像样本是有54个波段（这是经过特征提取的，每个波段代表一个特征信息），最后一个波段是标签，因此我的目标是利用前53个特征波段来预测最后一个标签波段。现在我有一组大的遥感影像作为训练样本，但是这些图像形状不规则，因此想要输入深度学习模型中进行训练就需要进行图像切割，变成比如64*64的小尺寸的多个样本，这个就是本文想要分享的内容。

        废话不多说，直接开始上码。

二、异常值处理

        首先导入要用到的python库

import os
from osgeo import gdal
import numpy as np
import pandas as pd

  这里我进行了一步异常值处理，由于我的第54波段，即标签存在一定的问题，因此我进行了一步根据阈值的滤除，去除了一些极大和极小值

 先定义输入和输出的文件夹路径（输出文件夹一定要先创建好，不然会报错），并获取这个路径下的tif文件，如果大家是单个tif文件，可以直接将tiff_files改为文件路径

# 输入文件夹路径，包含要处理的TIFF文件，即大遥感图像
input_folder = "UNet实验数据/原始数据"

# 输出文件夹路径，用于存储处理后的TIFF文件
output_folder = "UNet实验数据/异常值处理"

# 获取输入文件夹中的所有TIFF文件
tiff_files = [f for f in os.listdir(input_folder) if f.endswith('.tif')]

下面写一段循环来处理这个文件夹下的遥感图像异常值，if语句后是处理的条件，这里我写的channel == 53，即第54波段，如果存在异常值，即 <label_min 或 >label_max，就变为nan值，这里label_min和label_max是我之前定义的，大家自行可以换成别的值。

# 循环处理每个TIFF文件
for tiff_file in tiff_files:
    # 构建输入和输出文件的完整路径
    input_path = os.path.join(input_folder, tiff_file)
    output_path = os.path.join(output_folder, tiff_file)

    # 打开输入TIFF文件
    ds = gdal.Open(input_path, gdal.GA_ReadOnly)

    # 获取输入TIFF文件的宽度、高度和通道数
    image_width = ds.RasterXSize
    image_height = ds.RasterYSize
    num_channels = ds.RasterCount

    # 创建输出TIFF文件
    driver = gdal.GetDriverByName("GTiff")
    out_ds = driver.Create(output_path, image_width, image_height, num_channels, gdal.GDT_Float32)

    # 处理每个通道
    for channel in range(num_channels):
        channel_data = ds.GetRasterBand(channel + 1).ReadAsArray()

        if channel == 53:
            channel_data = np.where((channel_data < label_min) | (channel_data > label_max), np.nan, channel_data)
            
        # 写入处理后的通道数据
        out_ds.GetRasterBand(channel + 1).WriteArray(channel_data)

    # 释放数据集
    out_ds = None

print("处理完成。")

        到这里就完成了异常值的处理，在输出文件夹中会出现与原始数据相同的 处理完异常值后的影像。

三、空值处理

        由于我的样本中有一些波段的某些像元是nan，但是其他波段不是nan，因此需要对这些像元进行处理。比如，如果第10波段的第一行第一列像元值是nan，但是其他波段的第一行第一列像元都有值，那么就把所有波段第一行第一列的像元值都变为nan，即波段中的空值处理。具体代码如下（这里也需要注意的是，输出文件夹需要提前创建）：

# 输入文件夹路径，包含要处理的TIFF文件
input_folder = "UNet实验数据/异常值处理"
# 输出文件夹路径，用于存储处理后的TIFF文件
output_folder = "UNet实验数据/空值处理"

# 获取输入文件夹中的所有TIFF文件
tiff_files = [f for f in os.listdir(input_folder) if f.endswith('.tif')]

# 循环处理每个TIFF文件
for tiff_file in tiff_files:
    # 构建输入和输出文件的完整路径
    input_path = os.path.join(input_folder, tiff_file)
    output_path = os.path.join(output_folder, tiff_file)

    # 打开输入TIFF文件
    ds = gdal.Open(input_path, gdal.GA_ReadOnly)

    # 获取输入TIFF文件的宽度、高度和通道数
    image_width = ds.RasterXSize
    image_height = ds.RasterYSize
    num_channels = ds.RasterCount

    # 创建输出TIFF文件
    driver = gdal.GetDriverByName("GTiff")
    out_ds = driver.Create(output_path, image_width, image_height, num_channels, gdal.GDT_Float32)

    # 处理每个通道
    for channel in range(num_channels):
        channel_data = ds.GetRasterBand(channel + 1).ReadAsArray()

        # 创建一个掩码，用于检测NaN值
        nan_mask = np.isnan(channel_data)

        # 将NaN值应用到所有波段
        for other_channel in range(num_channels):
            other_channel_data = ds.GetRasterBand(other_channel + 1).ReadAsArray()
            other_channel_data[nan_mask] = np.nan
            out_ds.GetRasterBand(other_channel + 1).WriteArray(other_channel_data)

    # 释放数据集
    out_ds = None

print("处理完成。")

四、零值（均值）填充

     由于深度学习没法处理nan值，因此这里我单独用一段代码，将图像中的nan值用0值（或波段均值）来代替，思路和上面基本一样。

        在语句 channel_data[nan_mask] = 0 下面我还加了一句 channel_data[nan_mask] = np.nanmean(channel_data)，这是分别用0值和波段均值填充，大家可以根据情况更改。

# 输入文件夹路径，包含要处理的TIFF文件
input_folder = "UNet实验数据/空值处理"

# 输出文件夹路径，用于存储处理后的TIFF文件
output_folder = "UNet实验数据/零值填充"

# 获取输入文件夹中的所有TIFF文件
tiff_files = [f for f in os.listdir(input_folder) if f.endswith('.tif')]

# 循环处理每个TIFF文件
for tiff_file in tiff_files:
    # 构建输入和输出文件的完整路径
    input_path = os.path.join(input_folder, tiff_file)
    output_path = os.path.join(output_folder, tiff_file)

    # 打开输入TIFF文件
    ds = gdal.Open(input_path, gdal.GA_ReadOnly)

    # 获取输入TIFF文件的宽度、高度和通道数
    image_width = ds.RasterXSize
    image_height = ds.RasterYSize
    num_channels = ds.RasterCount

    # 创建输出TIFF文件
    driver = gdal.GetDriverByName("GTiff")
    out_ds = driver.Create(output_path, image_width, image_height, num_channels, gdal.GDT_Float32)

    # 处理每个通道
    for channel in range(num_channels):
        channel_data = ds.GetRasterBand(channel + 1).ReadAsArray()

        # 创建一个掩码，用于检测NaN值
        nan_mask = np.isnan(channel_data)

        # 用0值/均值填充NaN值
        channel_data[nan_mask] = 0
        # channel_data[nan_mask] = np.nanmean(channel_data)

        out_ds.GetRasterBand(channel + 1).WriteArray(channel_data)

    # 释放数据集
    out_ds = None

print("处理完成。")

五、滑动窗口分割

        下面就是滑窗分割了，分割思想很简单，本质也是为了增加样本量：

       如图所示，比如我要切割成32*32的图像，间隔（步长）16，首先在第一个像素位置（第一行第一列）进行分割，选取长度和宽度为32的窗口，切取该窗口内的图像所有像素，形成第一个样本；然后将这个32*32的窗口按照步长（假如为16）移动至第二个像素位置（第一行第十七列），在该位置切取该窗口内的图像所有像素，形成第二个样本；以此类推，完成整幅大的遥感影像的切割，实现样本制作。

直接上码：

这里需要改的就是路径，x,y的间隔（即步长），小尺寸子图像大小

# 输入文件夹路径，包含要分割的TIFF文件
input_folder = "UNet实验数据/零值填充"

# 输出文件夹路径，用于存储分割后的子图像
output_folder = "UNet实验数据/图像切割16"

# 定义x和y的间隔
x_interval = 16  # 横向间隔
y_interval = 16  # 纵向间隔

# 定义子图像的大小
sample_width = 32  # 设置子图像宽度为32像素
sample_height = 32  # 设置子图像高度为32像素

# 获取输入文件夹中的所有TIFF文件
tiff_files = [f for f in os.listdir(input_folder) if f.endswith('.tif')]

# 初始化计数器
sample_counter = 0

# 循环处理每个TIFF文件
for tiff_file in tiff_files:
    # 构建输入和输出文件的完整路径
    input_path = os.path.join(input_folder, tiff_file)
    
    # 打开输入TIFF文件
    ds = gdal.Open(input_path)
    
    # 获取输入TIFF文件的宽度、高度和通道数
    image_width = ds.RasterXSize
    image_height = ds.RasterYSize
    num_channels = ds.RasterCount

    # 循环分割图像
    for y in range(0, image_height - sample_height + 1, y_interval):
        for x in range(0, image_width - sample_width + 1, x_interval):
            # 切割子图像
            sample = np.zeros((num_channels, sample_height, sample_width), dtype=np.float32)
            for channel in range(num_channels):
                channel_data = ds.GetRasterBand(channel + 1).ReadAsArray(x, y, sample_width, sample_height)
                sample[channel, :, :] = channel_data

            # 检查子图像的尺寸是否符合要求
            if sample.shape == (num_channels, sample_height, sample_width):
                # 构建输出文件路径
                filename = os.path.join(output_folder, f'sample_{sample_counter}.tif')

                # 保存子图像为TIFF文件
                driver = gdal.GetDriverByName("GTiff")
                out_ds = driver.Create(filename, sample_width, sample_height, num_channels, gdal.GDT_Float32)
                for channel in range(num_channels):
                    out_ds.GetRasterBand(channel + 1).WriteArray(sample[channel, :, :])
                out_ds = None

                # 增加计数器
                sample_counter += 1


print("处理完成。")

最终就能生成滑窗分割后的图像了

六、数据增强

        有些朋友还需要进行数据增强，这里我就一起写了，但是要注意的是，在分割和增强之前一定要准备好测试集，不然容易信息重复，导致模型训练和测试精度过高，测试不准。

        这里需要更改的是输入和输出文件夹路径，目标图像尺寸和波段数量。

# 输入和输出文件夹路径
input_folder = "UNet实验数据/图像切割16"
output_folder = "UNet实验数据/数据增强16"

# 定义目标图像尺寸
target_height = 32
target_width = 32
num_bands =  54 # 波段数量

# 获取输入文件夹中的所有.tif文件
tiff_files = [file for file in os.listdir(input_folder) if file.endswith(".tif")]

# 循环处理每个TIFF文件
for tiff_file in tiff_files:
    # 构建输入文件的完整路径
    input_path = os.path.join(input_folder, tiff_file)
    
    # 打开TIFF文件
    ds = gdal.Open(input_path)

    if ds is None:
        raise Exception("无法打开TIFF图像")

    # 获取图像的宽度、高度
    image_width = ds.RasterXSize
    image_height = ds.RasterYSize

    # 读取图像数据
    image_data = ds.ReadAsArray()

    # 水平翻转
    flipped_horizontal = np.fliplr(image_data)
    output_file_name = f"{os.path.splitext(tiff_file)[0]}_flipped_horizontal.tif"
    output_path = os.path.join(output_folder, output_file_name)
    driver = gdal.GetDriverByName("GTiff")
    out_ds = driver.Create(output_path, target_width, target_height, num_bands, gdal.GDT_Float32)
    out_ds.SetProjection(ds.GetProjection())
    out_ds.SetGeoTransform(ds.GetGeoTransform())
    for band_idx in range(num_bands):
        out_ds.GetRasterBand(band_idx + 1).WriteArray(flipped_horizontal[band_idx])
    out_ds = None

    # 垂直翻转
    flipped_vertical = np.flipud(image_data)
    output_file_name = f"{os.path.splitext(tiff_file)[0]}_flipped_vertical.tif"
    output_path = os.path.join(output_folder, output_file_name)
    out_ds = driver.Create(output_path, target_width, target_height, num_bands, gdal.GDT_Float32)
    out_ds.SetProjection(ds.GetProjection())
    out_ds.SetGeoTransform(ds.GetGeoTransform())
    for band_idx in range(num_bands):
        out_ds.GetRasterBand(band_idx + 1).WriteArray(flipped_vertical[band_idx])
    out_ds = None

    # 保存原始图像
    output_file_name = f"{os.path.splitext(tiff_file)[0]}_original.tif"
    output_path = os.path.join(output_folder, output_file_name)
    out_ds = driver.Create(output_path, target_width, target_height, num_bands, gdal.GDT_Float32)
    out_ds.SetProjection(ds.GetProjection())
    out_ds.SetGeoTransform(ds.GetGeoTransform())
    for band_idx in range(num_bands):
        out_ds.GetRasterBand(band_idx + 1).WriteArray(image_data[band_idx])
    out_ds = None

print("数据增强完成")

七、效果展示

        这里给大家看一下原始图像和分割后的图像对比：

原始图像：

异常值处理，空值处理（左，Disp #1），均值填充（右，Disp #2）的图像

图像切割后的图像：

        至此，就可以用处理好的二维图像样本输入深度学习模型了，有时间我也会写一篇如何用这些样本进行遥感图像深度学习语义分割的完整过程，包括数据处理，特征分析，深度学习模型构建，训练，验证，精度评价，结果可视化等等。

# python # 深度学习 # 数据挖掘
Logo
腾讯云开发者社区

腾讯云面向开发者汇聚海量精品云计算使用和开发经验，营造开放的云计算技术生态圈。

更多推荐

JAVA使用 FtpClient连接频繁出现 java.net.SocketTimeoutException: Accept timed out

问题描述：我这边的系统通过FTP访问另一台服务器的目录进行文件读取和发送，但是总是频繁的出现 java.net.SocketTimeoutException: Accept timed out超时现象。但是该异常对系统正常运行并没有什么影响。只是在log日志中多了这些错误日志，如下：2020-10-14 04:03:15.028 ERROR pool-4-thread-2 FTPUtils:301

avatar 腾讯云开发者社区

python QT 图片缩放，移动

利用python的QT库搭建图片显示界面，实现图片的拖拽、以鼠标中心缩放等功能。

avatar 腾讯云开发者社区
cover

阿里云服务器搭建frps实现内网穿透

avatar 腾讯云开发者社区