碧
碧海苍梧
V1
2023/05/26阅读:8主题:默认主题
遥感影像深度学习样本对制作教程4——输入和标签空间范围不一致的情况
上一篇文章介绍了遥感影像深度学习样本对制作的常用技巧,但实际使用中并不是所有数据都是那么的规整的,比如输入数据和标签可能存在空间范围不一致的情况。关注公众号GeodataAnalysis
,回复20230526获取示例数据和代码,这个系列的代码都放在一起,上手运行一下代码更容易弄懂。
当输入数据和标签空间范围不一致时,需要考虑二者的像元在空间上能否一一对应,再分别采取不同的方式处理数据:
-
能一一对应时:直接提取二者的重叠范围,根据重叠部分制作标签。 -
不能一一对应时:重采样使其一一对应。
所以总的来说关键的步骤只有三点,提取重叠范围,根据重叠部分采用教程1的方法制作标签,以及重采样。下面分别介绍这三点,想具体了解可以自己下载代码和数据运行一下,这里只介绍关键步骤,具体的不多介绍。
1 计算两幅影像的重叠范围
这里需要的两幅影像必须是像元能够一一对应的,如示例数据中LC08_L1TP_018030_20150907_20200908_02_T1_B1.TIF
和LC08_L1TP_018030_20150923_20200908_02_T1_B1.TIF
就是空间范围不一致,但像元一一对应的示例。
def overlap(src1, src2):
gt = src1.transform
w, h = src1.width, src1.height
src1_extent = (gt[2], gt[2]+gt[0]*w, gt[5], gt[5]+gt[4]*h)
gt = src2.transform
w, h = src2.width, src2.height
src2_extent = (gt[2], gt[2]+gt[0]*w, gt[5], gt[5]+gt[4]*h)
mask_extent = (
max((src1_extent[0], src2_extent[0])), # west
min((src1_extent[1], src2_extent[1])), # east
min((src1_extent[2], src2_extent[2])), # north
max((src1_extent[3], src2_extent[3])) # south
)
min_row, min_col = src1.index(mask_extent[0], mask_extent[2])
max_row, max_col = src1.index(mask_extent[1], mask_extent[3])
src1_window = Window.from_slices((min_row, max_row), (min_col, max_col))
min_row, min_col = src2.index(mask_extent[0], mask_extent[2])
max_row, max_col = src2.index(mask_extent[1], mask_extent[3])
src2_window = Window.from_slices((min_row, max_row), (min_col, max_col))
src1_array = src1.read(1, window=src1_window)
src2_array = src2.read(1, window=src2_window)
src1_mask = np.ones(shape=src1_array.shape[-2:])
mask = src1_array==src1.nodata
src1_mask[mask] = 0
src2_mask = np.ones(shape=src2_array.shape[-2:])
mask = src2_array==src2.nodata
src2_mask[mask] = 0
mask = (src1_mask == 1) & (src2_mask == 1)
transform = from_bounds(mask_extent[0], mask_extent[-1], mask_extent[1],
mask_extent[-2], mask.shape[1], mask.shape[0])
return src1_window, src2_window, mask, transform
2 根据重叠范围制作样本
数据准备。
ids = ['LC08_L1TP_018030_20150907_20200908_02_T1', 'LC08_L1TP_018030_20150923_20200908_02_T1']
id_types = ['img', 'label']
img_src = rio.open(f'./data/{ids[0]}/{ids[0]}_B1.TIF')
label_src = rio.open(f'./data/{ids[1]}/{ids[1]}_B1.TIF')
img_window, label_window, mask, transform = overlap(img_src, label_src)
样本制作。
image_size, slide = 256, 256
# 用于记录每个影像块的位置
df = pd.DataFrame(columns=['patch', 'id', 'row', 'col', 'image_size', 'type'])
patch_num = 0
patch_pair = 1
height, width = mask.shape
for row in range((height - image_size) // slide + 1):
for col in range((width - image_size) // slide + 1):
for window, id, id_type in zip([img_window, label_window], ids, id_types):
patch_row = window.row_off + row*slide
patch_col = window.col_off + col*slide
mask_patch = mask[row*slide: row*slide + image_size,
col*slide: col*slide + image_size]
if not np.all(mask_patch):
break
df.loc[patch_num, 'patch'] = patch_pair
df.loc[patch_num, 'type'] = id_type
df.loc[patch_num, 'id'] = id
df.loc[patch_num, 'row'] = patch_row
df.loc[patch_num, 'col'] = patch_col
if patch_num % 2 == 1:
patch_pair += 1
patch_num += 1
df['image_size'] = image_size
df.head()
3 重采样
示例数据1分辨率为30米,示例数据2分辨率为480米(教程2中下载的土地覆盖数据),如二者像元是一一对应的话,示例的一个像元会覆盖示例1的16个像元,但目前数据存在偏差,所以先对示例2进行重采样。
读取数据。
src1 = rio.open(path1)
src2 = rio.open(path2)
获取示例1的坐标信息。
gt = src1.transform
w, h = src1.width, src1.height
src1_extent = (gt[2], gt[2]+gt[0]*w, gt[5], gt[5]+gt[4]*h)
计算示例二重采样后的行列数。
src2_res = src2.transform[0]
src2_w = (src1_extent[1] - src1_extent[0]) // src2_res
src2_h = (src1_extent[2] - src1_extent[3]) // src2_res
计算示例2重采样后的空间范围。
src2_extent = [src1_extent[0], src1_extent[0]+src2_w*src2_res,
src1_extent[2], src1_extent[2]-src2_h*src2_res]
计算示例2重采样后的GeoTransform,也就是仿射变换六参数。
dst_transform, dst_width, dst_height = calculate_default_transform(
src_crs=src2.crs,
dst_crs=src2.crs,
width=src2_w,
height=src2_h,
left=src2_extent[0],
bottom=src2_extent[3],
right=src2_extent[1],
top=src2_extent[2]
)
执行重采样,并保存为tif。
dst_array = np.empty((dst_height, dst_width), dtype=np.uint8)
reproject(
# 源文件参数
source=src2.read(1),
src_crs=src2.crs,
src_transform=src2.transform,
# 目标文件参数
destination=dst_array,
dst_transform=dst_transform,
dst_crs=src2.crs,
num_threads=1
)
array_to_tif('reproject.tif', dst_array, src2.crs, dst_transform)
作者介绍
碧
碧海苍梧
V1