I have an xarray dataset that I want to convert to image in for a WMS server I am creating in python. The problem is that my when I display my WMS with OpenLayer in my frontend, the tiles seems shifted.

I wonder if my reprojection and/or interpolation are the cause of this shift.
OpenLayers always request 256px by 256 px tiles .

My data have much more point that 256 in each axis, whatever the zoom level, so my interpolation is subsampling.

Her is my code:

def find_closest_indices(data, x_min, x_max, y_min, y_max):
    # Find the closest indices that are just outside the BBOX
    x_indices = [np.argmin(np.abs(data.x.values - x_min)),
                 np.argmin(np.abs(data.x.values - x_max))]
    y_indices = [np.argmin(np.abs(data.y.values - y_max)),
                 np.argmin(np.abs(data.y.values - y_min))]

    print(f"X indices: {x_indices}")
    print(f"Y indices: {y_indices}")

    # Extend indices to include the next point outside the BBOX
    x_indices = [max(0, x_indices[0] - 1),
                 min(len(data.x), x_indices[1] + 1)]
    y_indices = [max(0, y_indices[1] - 1),
                 min(len(data.y), y_indices[0] + 1)]

    return x_indices, y_indices


def extract_bbox(bbox_str: str) -> tuple:
    """Extract bounding box coordinates from a string"""
    return tuple(map(float, bbox_str.split(",")))


def transform_coordinates(proj_in, proj_out, lon_min, lat_min, lon_max, lat_max):
    """Transform bounding box coordinates to match dataset's coordinate system"""
    x_min, y_min = transform(proj_in, proj_out, lon_min, lat_min)
    x_max, y_max = transform(proj_in, proj_out, lon_max, lat_max)
    return x_min, y_min, x_max, y_max


def normalize_data(data: np.ndarray) -> np.ndarray:
    """Normalize data values to fit within a fixed range"""
    return (data - COLOR_SCALE_MIN) / (COLOR_SCALE_MAX - COLOR_SCALE_MIN)


# @lru_cache(maxsize=128)
def get_interpolated_data(filepath, var_name, bbox, width, height, crs, padding=1):
    """Extract and interpolate data from a NetCDF file"""
    lon_min, lat_min, lon_max, lat_max = bbox
    # Define projections (example UTM zone 33, adjust as needed)
    x_min, y_min, x_max, y_max = transform_coordinates(
        crs, "EPSG:32631", *bbox)

    # print(f"X min: {x_min}, X max: {x_max}")
    # print(f"Y min: {y_min}, Y max: {y_max}")

    ds = xr.open_dataset(filepath)
    data = ds[var_name]

    x_indices, y_indices = find_closest_indices(
        data, x_min, x_max, y_min, y_max)

    print(f"Data shape: {data.x.shape}, {data.y.shape}")

    # Select the first wavelength layer
    subset = data.isel(wl=2).isel(
        x=slice(x_indices[0], x_indices[1]),
        y=slice(y_indices[1], y_indices[0])
    )
    print(
        f"Subset shape: {subset.shape}, x: {subset.x.shape}, y: {subset.y.shape}")

    data_np = subset.values
    print(f"Data shape: {data_np.shape}")

    # Create a mask for NaN values
    nan_mask = np.isnan(data_np)
    print(f"NaN mask: {nan_mask.sum()} NaNs")

    # Apply NaN mask to keep NaNs in the data

    if data_np.shape == (0, 0) or data_np.shape[0] == 0 or data_np.shape[1] == 0:
        return data_np

    data_normalized = cv2.normalize(
        data_np, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)
    data_normalized[nan_mask] = 0
    # Normalize data to a range suitable for OpenCV (e.g., uint8)
    # data_lon_min, data_lat_min, data_lon_max, data_lat_max = transform_coordinates(
    #     "EPSG:32631", crs, subset.x.min().item(), subset.y.min().item(), subset.x.max().item(), subset.y.max().item())

    # Convert bbox and data extent to pixel coordinates
    lon_range = x_max - x_min
    lat_range = x_max - x_min

    data_x_min = subset.x.min().item()
    data_y_min = subset.y.min().item()
    data_x_max = subset.x.max().item()
    data_y_max = subset.y.max().item()

    print(f"Data lon min: {data_x_min}, lon max: {data_x_max}")
    print(f"Data lat min: {data_y_min}, lat max: {data_y_max}")

    print(f"lon min: {x_min}, lon max: {x_max}")
    print(f"lat min: {y_min}, lat max: {y_max}")

    print((f"Range: {lon_range}, {lat_range}"))

    # Calculate the pixel coordinates of the data
    data_px_min = (data_x_min - x_min) / (x_max - x_min) * width
    data_py_min = (y_max - data_y_max) / (y_max - y_min) * height
    data_px_max = (data_x_max - x_min) / (x_max - x_min) * width
    data_py_max = (y_max - data_y_min) / (y_max - y_min) * height

    print("Before rounding")
    print(f"Data x min: {data_px_min}, x max: {data_px_max}")
    print(f"Data y min: {data_py_min}, y max: {data_py_max}")

    # Parse the pixel coordinates to integers
    data_px_min = int(data_px_min)
    data_py_min = int(data_py_min)
    data_px_max = int(data_px_max) + 1
    data_py_max = int(data_py_max) + 1

    if data_px_max > width:
        data_px_max = width
    if data_py_max > height:
        data_py_max = height
    if data_px_min < 0:
        data_px_min = 0
    if data_py_min < 0:
        data_py_min = 0

    print("After rounding")
    print(f"Data x min: {data_px_min}, x max: {data_px_max}")
    print(f"Data y min: {data_py_min}, y max: {data_py_max}")

    print(f"x_min: {x_min}, x_max: {x_max}")
    print(f"y_min: {y_min}, y_max: {y_max}")

    # Calculate the width and height of the data in pixel coordinates
    data_pixel_width = data_px_max - data_px_min
    data_pixel_height = data_py_max - data_py_min

    print(f"Data pixel width: {data_pixel_width}, height: {data_pixel_height}")

    data_zeros = np.zeros(
        (data_pixel_height, data_pixel_width), dtype=data_normalized.dtype)
    data

    # Resize the data to fit within its pixel coordinates
    data_resized = cv2.resize(
        data_normalized, (data_pixel_width, data_pixel_height), interpolation=cv2.INTER_AREA)

    print(f"Resized data shape: {data_resized.shape}")

    # Place the resized data into the target-sized array
    data_interpolated = np.zeros((height, width), dtype=data_resized.dtype)
    data_interpolated[data_py_min:data_py_max,
                      data_px_min:data_px_max] = data_resized

    data_interpolated = data_interpolated.astype(
        float)  # Convert to float to support NaN values
    data_interpolated[data_interpolated == 0] = np.nan

    return data_interpolated

def create_color_mapped_image(data, width, height, cmap_name='plasma'):
    """Create a color-mapped image from normalized data"""
    if data.size == 0:
        print("Data is empty")
        # Return a transparent image of default size
        return Image.new('RGBA', (256, 256), (255, 255, 255, 0))
    norm = Normalize(vmin=0, vmax=1)
    # print(f"Normalized data min: {norm.vmin}, max: {norm.vmax}")
    cmap = plt.get_cmap(cmap_name)
    mappable = ScalarMappable(norm=norm, cmap=cmap)
    mapped_data = (mappable.to_rgba(data, bytes=True)).astype(np.uint8)

    image = Image.fromarray(mapped_data, mode="RGBA")
    return image

Does anyone as experience in interpolation for tiles images ?