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