|
| 1 | +import dask.array |
| 2 | +import numpy as np |
| 3 | + |
| 4 | +import xarray |
| 5 | + |
| 6 | +x = np.arange(56) |
| 7 | +y = np.arange(100, 119) |
| 8 | +t = np.arange(2) |
| 9 | +z = np.arange(9) |
| 10 | + |
| 11 | +# The first array |
| 12 | +projected = xarray.DataArray( |
| 13 | + dask.array.full((x.size, y.size), np.nan, chunks=(10, 10)), |
| 14 | + coords={"x": x, "y": y}, |
| 15 | + name="projected", |
| 16 | +) |
| 17 | + |
| 18 | +# The array whose values will be used |
| 19 | +original = xarray.DataArray( |
| 20 | + np.full((t.size, z.size), 12.345), coords={"t": t, "z": z}, dims=["t", "z"] |
| 21 | +).chunk({"t": 10}) |
| 22 | + |
| 23 | +t2x = np.linspace(0, x.size, t.size) |
| 24 | +z2y = np.linspace(y.min(), y.max(), z.size) |
| 25 | +tmp = t2x[..., None] @ z2y[None, ...] |
| 26 | + |
| 27 | +# The array with the coordinate correspondence |
| 28 | +pp = np.full((t.size, z.size, 2), np.nan) |
| 29 | +pp[..., 0] = tmp / z2y.max() |
| 30 | +pp[..., 1] = tmp / t2x.max() |
| 31 | + |
| 32 | +pixel_positions = xarray.DataArray( |
| 33 | + pp, coords={"t": t, "z": z, "nav": ["x", "y"]}, dims=["t", "z", "nav"] |
| 34 | +).chunk({"t": 10}) |
| 35 | + |
| 36 | + |
| 37 | +def match_coordinates(projected, pixel_positions): |
| 38 | + coordinates = projected.sel( |
| 39 | + x=pixel_positions[..., 0], y=pixel_positions[..., 1], method="nearest" |
| 40 | + ).coords |
| 41 | + return {"x": coordinates["x"], "y": coordinates["y"]} |
| 42 | + |
| 43 | + |
| 44 | +def place_pixel_on_grid(projected, pixel_positions, original): |
| 45 | + coordinates = match_coordinates(projected, pixel_positions) |
| 46 | + print(f"{coordinates=}") |
| 47 | + print(f"{projected.shape} -- {original.shape}") |
| 48 | + print(f"pixel_positions coords {pixel_positions.coords}") |
| 49 | + projected_chunk = projected.copy() |
| 50 | + projected_chunk.loc[{"x": coordinates["x"], "y": coordinates["y"]}] = original |
| 51 | + return projected_chunk |
| 52 | + |
| 53 | + |
| 54 | +xarray.map_blocks( |
| 55 | + place_pixel_on_grid, projected, args=[pixel_positions, original], template=projected |
| 56 | +).compute() |
0 commit comments