Clip grid

import geopandas as gp
import xarray as xr

from dask_jobqueue import PBSCluster
from dask.distributed import Client

cluster = PBSCluster(
    cores=80,
    processes=40,
    memory='185GB',
    interface='ib0',
    local_directory='$TMPDIR',
    queue='development',
    walltime='06:00:00',
    job_script_prologue=["spack env activate analysis"]
)
print(cluster.job_script())

cluster.scale(jobs=25)
client = Client(cluster)
client

# these are CHM ugrid regridded with chm_ugrid2grid
all = []
for y in range(2016, 2024):
    ds_y = xr.open_zarr(
        f"../monthly_mean_{y}.zarr",
        chunks="auto",
        consolidated=False,
        decode_cf=True,
    )
    all.append(ds_y)

# ensure alifgn ment
aligned = xr.align(*all, join="override", exclude="time")
ds = xr.combine_by_coords(aligned, combine_attrs="override")

# ensure spatial dims are clearly marked for rioxarray
ds = ds.rio.set_spatial_dims(x_dim="longitude", y_dim="latitude", inplace=False)
ds = ds.rio.write_crs("EPSG:4326", inplace=False)

# clip
shp = gp.read_file("myclipping.shp")
shp = shp.to_crs("epsg:4326")

clipped = ds.rio.clip(shp.geometry, drop=True)

# remove weird Zarr/GeoZarr encodings that confuse NetCDF writers
clipped.encoding = {}
for name, var in clipped.variables.items():
    print( var.attrs,var.encoding )
    var.attrs = {}
    var.encoding = {}

# dim order
clipped = clipped.transpose("time", "latitude", "longitude")

clipped.to_netcdf(
    "monthly_2016-2024_swe_sd.nc",
    engine="netcdf4",   # or omit engine= and let xarray choose
)