遥感影像深度学习样本对制作教程1

本文结合实例详细讲解了如何使用Python制作遥感影像深度学习样本对，关注公众号GeodataAnalysis，回复20230427获取示例数据和代码。

受到计算机内存限制，深度学习算法无法直接对大幅影像和标签图像进行训练，因此需要对影像和标签图像进行再处理，将制作好的城市绿地影像和对应标签图像裁剪成一系列相同尺寸的小图像放入模型训练。

常用的裁剪方法有规则格网裁剪和滑动窗口裁剪。规则格网裁剪，也称棋盘裁剪，顾名思义就是将影像裁剪成棋盘形状固定尺寸的小图像，该方法得到的图像比较规则，但是缺点是获取的图像数量较少，如下图a所示；滑动窗口裁剪在格网裁剪的基础上加以改进，通过设置指定的滑动间距进行滑动裁剪。滑动窗口裁剪相比规则格网裁剪可以获得更多的样本数量，能够充分利用城市绿地样本数据集，如下图b所示。

可见，只要将滑动窗口裁剪的滑动间隔设为和图像尺寸相同，那么滑动窗口裁剪就等同于规则裁剪。因此，本文只介绍滑动窗口裁剪。

首先导入必须的包，其中pandas用于存储每个裁剪影像块（patch）的信息，rasterio用于读取和保存栅格影像数据。

import os
import numpy as np
import pandas as pd
import rasterio as rio
from rasterio.windows import Window
from rasterio.transform import from_bounds

确定裁剪影像块的信息

遥感影像深度学习用到的数据一般是公开的，为了方便把自己训练模型的所用的数据分享给别人，不用分享整个数据集，记录每个影像块的信息，分享这个信息就很方便。而且，这也方便后续处理，如从这些影像块中进行筛选去除不合格的影像块。

首先确定影像块尺寸，和滑动间隔。

image_size = 256
slide = 256

读取输入数据和样本数据，判断二者的transform是否相同。若相同，则说明其空间范围一直，分辨率一致，最后记录二者的宽和高。测试数据为landsat 8数据，各个波段分开存储，输入数据和标签数据用的是同一个，仅作为示例。id相当于这副影像的身份信息，可以利用它通过公开渠道下载这幅影像。

# 打开模型输入数据
id1 = 'LC08_L1TP_018030_20150907_20200908_02_T1'
src1 = rio.open(f'./data/{id1}/{id1}_B1.TIF')
# 打开标签数据
id2 = 'LC08_L1TP_018030_20150907_20200908_02_T1'
src2 = rio.open(f'./data/{id2}/{id2}_B1.TIF')

assert src1.transform == src2.transform
# 读取影像的高度和宽度
height, width = src1.height, src1.width

输入数据往往具有nodata值，在训练模型我们有时并不想要输入数据包含nodata值，这时候就需要先创建一个掩膜，记录哪些像元是nodata，用于后续的裁剪。

# 根据输入数据创建掩膜
mask = np.full(shape=(height, width), fill_value=False)
data = src1.read(1)
# 所有空值区域设为True
mask[data==src1.nodata] = True

创建一个pandas的DataFrame，用于记录影像块的信息。

df = pd.DataFrame(columns=['patch', 'id', 'row', 'col', 'image_size', 'type'])

计算每个影像块在整幅影像的相对位置，即影像块的左上角的像元位于整幅影像的第几行第几列。并判断这个影像块是否包含nodata值，决定是否保留。

patct_num = 1
num = 0

for row in range((height - image_size) // slide + 1):
    for col in range((width - image_size) // slide + 1):
        # 读取每个patch的在输入数据的nodata掩膜
        mask_patch = mask[row*slide: row*slide + image_size, 
                          col*slide: col*slide + image_size]
        # 若这个影像块包含nodata，舍弃
        if np.any(mask_patch):
            continue
        # 保存输入影像块
        df.loc[num, 'patch'] = patct_num
        df.loc[num, 'id'] = id1
        df.loc[num, 'row'] = row * slide
        df.loc[num, 'col'] = col * slide
        df.loc[num, 'image_size'] = image_size
        df.loc[num, 'type'] = 'img'
        num += 1
        # 保存标签影像块
        df.loc[num, 'patch'] = patct_num
        df.loc[num, 'id'] = id2
        df.loc[num, 'row'] = row * slide
        df.loc[num, 'col'] = col * slide
        df.loc[num, 'image_size'] = image_size
        df.loc[num, 'type'] = 'label'
        num += 1
        patct_num += 1

划分训练集、验证集、测试集

以下代码按照6：2：2的比例将数据划分为训练集、验证集和测试集。

train_ratio = 6
valid_ratio = 2
test_ratio = 2
split = ['train']*train_ratio*2 + ['val'] * valid_ratio*2 + ['test']*test_ratio*2
split = split * (df.shape[0] // (train_ratio + valid_ratio + test_ratio) + 1)
split = split[:df.shape[0]]
df['split'] = split

根据记录信息裁剪影像块

用于将numpy数组转换为GeoTiff文件保存的函数.

# numpy数组转tif
def array_to_tif(out_path, arr, crs, transform):
    # 获取数组的形状
    if arr.ndim==2:
        count = 1
        height, width = arr.shape
    elif arr.ndim==3:
        count = arr.shape[0]
        _, height, width = arr.shape
    else:
        raise ValueError
    
    with rio.open(out_path, 'w', 
                  driver='GTiff', 
                  height=height, width=width, 
                  count=count, 
                  dtype=arr.dtype, 
                  crs=crs, 
                  transform=transform) as dst:
        # 写入数据到输出文件
        if count==1:
            dst.write(arr, 1)
        else:
            for i in range(count):
                dst.write(arr[i, ...], i+1)

确定需要处理哪一个波段，目前测试数据的各个波段是分开的，后续的文章会介绍分波段存储的文件合并为一个文件。

band = 1

打开记录中的整幅影像，虽然测试中只用了一景影像，即输入和标签是相同的，但实际工作中可能会有很多景，相当于很多个id。注意，这里打开数据并不是将影像的数据读入内存，只是提供一个读取这幅影像数据的入口。

id_l8_srcs = {}
data_dir = './data'
for id in df['id'].unique():
    id_path = os.path.join(data_dir, id, id+f'_B{band}.TIF')
    id_l8_srcs[id] = rio.open(id_path)

创建训练集、验证集、测试集的输出文件夹，并在其中创建标签和输入数据的文件夹。

out_dir = './dataset'

for split in df['split'].unique():
    for type in df['type'].unique():
        split_dir = os.path.join(out_dir, split, type)
        os.makedirs(split_dir, exist_ok=True)

裁剪影像块，并将其保存到它该带待的位置。注意，这里保存的每个影像块都是具有空间坐标信息的，而不只是一个图片，这样就为一些预处理等操作提供了便利。

for patch in df.itertuples():
    l8_src = id_l8_srcs[patch.id]
    # 影像块的窗口
    win = Window.from_slices((patch.row, patch.row+patch.image_size), 
                             (patch.col, patch.col+patch.image_size))
    # 读取影像块的值
    l8_patch_array = l8_src.read(1, window=win)
    # 计算影像块的空间范围
    west, north = l8_src.xy(patch.row, patch.col, offset='ul')
    east, south = l8_src.xy(patch.row+patch.image_size, 
                            patch.col+patch.image_size, offset='ul')
    transform = from_bounds(west, south, east, north, 
                            patch.image_size, patch.image_size)
    # 保存影像块
    out_path = os.path.join(out_dir, patch.split, 
                            patch.type, f'{patch.patch}.TIF')
    array_to_tif(out_path, l8_patch_array, l8_src.crs, transform)

裁剪后的影像块示例：