迭代地按块加载图像，其中块部分重叠，这些块逐步加载。

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

然而，对于任何需要内核操作的处理（例如任何卷积滤波器，如平滑），这会导致在块边缘附近处理时出现问题。为了解决这个问题，每个块应该略微重叠，如下所示：

我的目标是按照以下方式加载块，其中重叠的数量可以根据需要进行调整。

到目前为止，我没有找到任何直接的方法来实现这一点。如有帮助，我将不胜感激。

- tachyon

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)

输出窗口：

- user2856

---

由于GDAL和rasterio都没有提供方便的方法来实现它，我想这是目前的解决方法。我在rasterio主组中开了一个帖子，讨论将此功能实现到rasterio中的问题。 - tachyon

@tachyon - 请查看错误修复编辑。原始版本从右侧和底部裁剪了一些像素。 - user2856