jgrss
commented
3 months ago Hey, can you post the input data so I can better see what you're trying to do?
Closed mmann1123 closed 3 months ago
jgrss
commented
3 months ago Hey, can you post the input data so I can better see what you're trying to do?
mmann1123
commented
3 months ago @jgrss Here's an example:
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)I think you can do the first approach if you set the overall bounds.
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)As for the VRT approach, that looks like it could work. However, there is a bug when passing the default integer EPSG -- gdal.BuildVRTOptions outputSRS does not like EPSG integers, so you could try setting the DataArray CRS attribute as a WKT string like src1.attrs["crs"] = src1.gw.crs_to_pyproj.to_wkt().
jgrss
commented
3 months ago Follow-up on the first option. In the case of more than two arrays, you might want to use numpy.
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)ok thanks I will play around with it. Looks like the VRT option is still throwing the error.
@jgrss Just wondering if you have some advice for stacking and mosaicing. I tried both of the following and am getting errors. Is there a way to do this?
error for stack & mosaic
Error for vrt: