0

I want to add a bunch of different raster maps in datasets with common data variables.

So far, I have been using my own complicated function to fix this problem (See below). But I have a suspicion this is not optimal, and I seem to be getting some weird artifacts, such as pixels being dragged out in lines. In the cases where no data overlaps, I understand I should be using the dxr1.merge(dxr2.interp_like(dxr2)) function, but this does not work in the case where I want to add the values in the data. Any suggestions? :)


def merge_raster(
    ds1: Union[xr.Dataset, xr.DataArray],
    ds2: Union[xr.Dataset, xr.DataArray],
    dataset: str = "area",
) -> xr.Dataset:
    """Merge two xarray raster-maps

    Args:
        ds1 (xr.Dataset): Dataset 1
        ds2 (xr.Dataset): Dataset 2
        dataset (str, optional): Name of data-variable to be added. Defaults to "area".

    Returns:
        xr.Dataset: Merged and added dataset
    """
    x1 = ds1.x.values
    x2 = ds2.x.values
    y1 = ds1.y.values
    y2 = ds2.y.values
    x = np.sort(np.unique(np.append(x1, x2)))
    y = np.sort(np.unique(np.append(y1, y2)))
    data = np.zeros((len(y), len(x)))

    ds1_data = ds1[dataset].values if isinstance(ds1, xr.Dataset) else ds1.values
    ds2_data = ds2[dataset].values if isinstance(ds2, xr.Dataset) else ds2.values

    xidx1 = x.searchsorted(x1)
    yidx1 = y.searchsorted(y1)
    data[yidx1[:, None], xidx1] += np.nan_to_num(ds1_data)

    xidx2 = x.searchsorted(x2)
    yidx2 = y.searchsorted(y2)
    data[yidx2[:, None], xidx2] += np.nan_to_num(ds2_data)

    data[data == 0] = np.nan
    data = xr.Dataset(
        data_vars={
            dataset: (
                [
                    "y",
                    "x",
                ],
                data,
            )
        },
        coords=dict(
            x=(
                ["x"],
                x,
            ),
            y=(
                ["y"],
                y,
            ),
        ),
    )
    return data

0 Answers0