在块部分重叠的情况下逐块迭代地加载图像

Question

试图处理大型卫星图像(~10GB).对于存储器有效处理,在每次迭代中将图像块(block/ tile)加载到存储器中.

示例代码如下:

def process_image(src_img, dst_img, band_id=1):
    with rasterio.open(src_img) as src:
        kwargs = src.meta
        tiles = src.block_windows(band_id)
        with rasterio.open(dst_img, 'w', **kwargs) as dst:
            for idx, window in tiles:
                print("Processing Block: ", idx[0]+1, ", ", idx[1]+1)
                src_data = src.read(band_id, window=window)
                dst_data = src_data ** 2 # Do the Processing Here
                dst.write_band( band_id, dst_data, window=window)
    return 0

然而,对于需要内核方式操作的任何类型的处理(例如任何convolve类似的滤波器smoothing),这导致在块的边缘附近的处理中的问题.为解决此问题,每个块应略微重叠,如下所示:

我的目标是按以下方式加载块,其中可以根据需要调整重叠量.

到目前为止,我没有找到任何直接的方法来实现这一目标.在这方面我会感激不尽.

Answer 1

你可以扩展你的窗口:

import rasterio as rio
from rasterio import windows

def overlapping_blocks(src, overlap=0, band=1):
    nols, nrows = src.meta['width'], src.meta['height']
    big_window = windows.Window(col_off=0, row_off=0, width=nols, height=nrows)
    for ji, window in src.block_windows(band):

        if overlap == 0:
            yield ji, window

        else:
            col_off = window.col_off - overlap
            row_off = window.row_off - overlap
            width = window.width + overlap * 2
            height = window.height + overlap * 2
            yield ji, windows.Window(col_off, row_off, width, height).intersection(big_window)

你会像这样使用它:

def process_image(src_img, dst_img, band_id=1):
    with rasterio.open(src_img) as src:
        kwargs = src.meta
        with rasterio.open(dst_img, 'w', **kwargs) as dst:
            for idx, window in overlapping_block_windows(src, overlap=1, band=band_id):
                print("Processing Block: ", idx[0]+1, ", ", idx[1]+1)
                src_data = src.read(band_id, window=window)
                dst_data = src_data ** 2 # Do the Processing Here
                dst.write_band( band_id, dst_data, window=window)
    return 0

这里有一种方法可以重叠块窗口和非块窗口,还有一个额外的参数boundless来指定是否将窗口扩展到栅格范围之外,对于无限读取非常有用:

from itertools import product
import rasterio as rio
from rasterio import windows


def overlapping_windows(src, overlap, width, height, boundless=False):
    """"width & height not including overlap i.e requesting a 256x256 window with 
        1px overlap will return a 258x258 window (for non edge windows)"""
    offsets = product(range(0, src.meta['width'], width), range(0, src.meta['height'], height))
    big_window = windows.Window(col_off=0, row_off=0, width=src.meta['width'], height=src.meta['height'])
    for col_off, row_off in offsets:

        window = windows.Window(
            col_off=col_off - overlap,
            row_off=row_off - overlap,
            width=width + overlap * 2,
            height=height + overlap * 2)

        if boundless:
            yield window
        else:
            yield window.intersection(big_window)


def overlapping_blocks(src, overlap=0, band=1, boundless=False):

    big_window = windows.Window(col_off=0, row_off=0, width=src.meta['width'], height=src.meta['height'])
    for ji, window in src.block_windows(band):

        if overlap == 0:
            yield window

        else:
            window = windows.Window(
                col_off=window.col_off - overlap,
                row_off=window.row_off - overlap,
                width=window.width + overlap * 2,
                height=window.height + overlap * 2)
            if boundless:
                yield window
            else:
                yield window.intersection(big_window)

输出窗口: