Stack & Mosaic Advice? about geowombat HOT 6 CLOSED

from geowombat.data import (
    l8_224077_20200518_B2,
    l8_224078_20200518_B2,
    l8_224077_20200518_B3,
    l8_224078_20200518_B3,
)

with gw.open([l8_224077_20200518_B2, l8_224077_20200518_B3], stack_dim="band") as src1:
    with gw.open(
        [l8_224078_20200518_B2, l8_224078_20200518_B3], stack_dim="band"
    ) as src2:

        with gw.open([src1, src2], mosaic=True) as src:
            print(src)

import dask.array as da
from geowombat.backends.rasterio_ import get_file_bounds

# Get the union of all bounding boxes
# Note that you probably want to pass `crs` in this method.
bounds = get_file_bounds([l8_224077_20200518_B2, l8_224078_20200518_B2], bounds_by="union", return_bounds=True)

with gw.config.update(ref_bounds=bounds):
    with gw.open([l8_224077_20200518_B2, l8_224077_20200518_B3], band_names=[1, 2], stack_dim="band") as src1:
        with gw.open([l8_224078_20200518_B2, l8_224078_20200518_B3], band_names=[1, 2], stack_dim="band") as src2:
            darray = da.maximum(src1, src2)
            print(darray)

import numpy as np

with gw.config.update(ref_bounds=bounds):
    with gw.open([l8_224077_20200518_B2, l8_224077_20200518_B3], band_names=[1, 2], stack_dim="band") as src1:
        with gw.open([l8_224078_20200518_B2, l8_224078_20200518_B3], band_names=[1, 2], stack_dim="band") as src2:
            darray = np.maximum.reduce(src1, src2, src3, src4, ...)
            print(darray)