Joining arrays, where boundary is not along an axis

[harryC-space-intelligence]

I am looking for a solution to mosaic arrays - by this I mean if I start with three arrays a, b, c, every value in my result will be taken either from a, b, or c.

The starting arrays overlap each other, and I want to define a boundary between where the result chooses one array over the other inside this region, based on some geometries that are not aligned to the x-y coordinates of the arrays.

Problem Illustration

1. My area of interest (AOI) is the green shaded region. The original data comes from regions indicated by the black outlines (satellite image footprints).

2. I generate three dataArrays that cover the overlapping regions between the AOI and the image footprints. They are all placed on a regular grid.

3. I want to join them such that the boundaries keep well away from the edges of the image footprints, eg:

It is not important exactly where the boundary lies, as long as every position in the final array is given a value from one of the starting arrays.

[scottyhq]

@harryC-space-intelligence thanks for the nice illustrations. If you're lucky enough to have satellite images with STAC metadata, the libraries stackstac and odc.stac have built-in mosaicing functions that take care of lots of details like ensuring your images are in the same CRS as well as using algorithms that work well with dask. If you don't have STAC metadata you might consider making it in order to work with these tools - it's pretty easy to create with rio-stac .

That said, if you already have a dataset you can create a mosaic using Xarray.reduce with something like the following:

def first_valid_value(arr, axis=None):
    idx = np.nanargmax(np.isfinite(arr), axis=axis, keepdims=True)
    vals = np.take_along_axis(arr, idx, axis=axis).squeeze() 
    return vals

mosaic = ds.value.reduce(first_valid_value, dim='landsat:scene_id')

I uploaded some code that walks through this with public data, and am also including just a couple plots below that show what this reduction does https://gist.github.com/scottyhq/ed8247f3ae1d42543f7bbfb02a5fa8ad

[harryC-space-intelligence]

Thanks @scottyhq for your answer.

Unfortunately my problem is slightly different (i.e. point 3: I want to join them such that the boundaries keep well away from the edges of the image footprints). If you look at my final illustration you can see how the join is equidistant from the edges of the images.

I found a way to implement this by looking for the indexes of null values along the rows of the data, something like this:

  def mask_along_overlap_center(left, right, dim='x'):

      left_finishes = left.bfill(dim=dim).isnull().argmax(dim=dim)
      right_starts = right.notnull().argmax(dim=dim)
      bisection = (left_finishes + right_starts) // 2 
  
      if left.chunks:
        ones = dask.array.ones(left.shape)
      else:
        ones = np.ones(left.shape)
  
      left_mask = xr.DataArray(ones, dims = left.dims, coords = left.coords).cumsum(dim=dim) < bisection
  
      return left_mask

[scottyhq]

Sure thing, always multiple ways to go about something. Thanks for sharing your working solution!