Skip to content

Bug: format_lon fails on latest xarray due to read-only coord values #65

@shenyulu

Description

@shenyulu

This is the error I discovered in the latest installed version. It seems that the lon array of xarray is immutable, and the following error occurred.

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[5], line 1
----> 1 data_u.regrid.linear(target_dataset)

File ~\anaconda3\envs\my\Lib\site-packages\xarray_regrid\regrid.py:44, in Regridder.linear(self, ds_target_grid, time_dim)
     33 """Regrid to the coords of the target dataset with linear interpolation.
     34 
     35 Args:
   (...)     41     Data regridded to the target dataset coordinates.
     42 """
     43 ds_target_grid = validate_input(self._obj, ds_target_grid, time_dim)
---> 44 ds_formatted = format_for_regrid(self._obj, ds_target_grid)
     45 return interp.interp_regrid(ds_formatted, ds_target_grid, "linear")

File ~\anaconda3\envs\my\Lib\site-packages\xarray_regrid\utils.py:276, in format_for_regrid(obj, target, stats)
    274 result = ensure_monotonic(result, coord)
    275 target = ensure_monotonic(target, coord)
--> 276 result = coord_handlers[coord_type]["func"](result, target, formatted_coords)
    278 # Coerce back to a single chunk if that's what was passed
    279 if isinstance(obj, xr.DataArray) and len(obj.chunksizes.get(coord, ())) == 1:

File ~\anaconda3\envs\my\Lib\site-packages\xarray_regrid\utils.py:387, in format_lon(obj, target, formatted_coords)
    386 if left_pad:
--> 387     lon_vals[:left_pad] = source_lon.values[-left_pad:] - 360
    388 if right_pad:
    389     lon_vals[-right_pad:] = source_lon.values[:right_pad] + 360

ValueError: assignment destination is read-only

The solution is quite simple, just add a copy() in utils.py to solve it.

def format_lon(
    obj: xr.DataArray | xr.Dataset, target: xr.Dataset, formatted_coords: dict[str, str]
) -> xr.DataArray | xr.Dataset:
    """Format the longitude coordinate by shifting the source grid to line up with
    the target anywhere in the range of -360 to 360, and then add a single wraparound
    padding column if the domain is inferred to be global and the east or west edges
    of the target lie outside the source grid centers.

    For example, with a source grid ranging from 0.5 to 359.5 and a target grid ranging
    from -180 to 180, the source grid would be shifted to -179.5 to 179.5 and then
    padded on both the left and right with wraparound values at -180.5 and 180.5 to
    provide full coverage for the target edge cells at -180 and 180.
    """
    lon_coord = formatted_coords["lon"]

    # Find a wrap point outside of the left and right bounds of the target
    # This ensures we have coverage on the target and handles global > regional
    source_vals = obj.coords[lon_coord].values
    target_vals = target.coords[lon_coord].values
    wrap_point = (target_vals[-1] + target_vals[0] + 360) / 2
    source_vals = np.where(
        source_vals < wrap_point - 360, source_vals + 360, source_vals
    )
    source_vals = np.where(source_vals > wrap_point, source_vals - 360, source_vals)
    obj = update_coord(obj, lon_coord, source_vals)

    obj = ensure_monotonic(obj, lon_coord)

    # Only pad if domain is global in lon
    source_lon = obj.coords[lon_coord]
    target_lon = target.coords[lon_coord]
    dx_s = source_lon.diff(lon_coord).max().values.item()
    dx_t = target_lon.diff(lon_coord).max().values.item()
    is_global_lon = source_lon.max().values - source_lon.min().values >= 360 - dx_s

    if is_global_lon:
        left_pad = (source_lon.values[0] - target_lon.values[0] + dx_t / 2) / dx_s
        right_pad = (target_lon.values[-1] - source_lon.values[-1] + dx_t / 2) / dx_s
        left_pad = int(np.ceil(np.max([left_pad, 0])))
        right_pad = int(np.ceil(np.max([right_pad, 0])))
        obj = obj.pad({lon_coord: (left_pad, right_pad)}, mode="wrap", keep_attrs=True)
        lon_vals = obj.coords[lon_coord].values.copy()    # <- add here!!!
        if left_pad:
            lon_vals[:left_pad] = source_lon.values[-left_pad:] - 360
        if right_pad:
            lon_vals[-right_pad:] = source_lon.values[:right_pad] + 360
        obj = update_coord(obj, lon_coord, lon_vals)
        obj = ensure_monotonic(obj, lon_coord)

    return obj

It would be great if there was time to fix this small mistake.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions