That would be super cool! Would need to decide if the flag_values are kept as well. I assume you would add this in da.cf.groupby(groups)?

A short example (using regionmask - could of course be simplified to only using xarray):

import cf_xarray
import regionmask
import xarray as xr

air = xr.tutorial.open_dataset("air_temperature")
mask = regionmask.defined_regions.ar6.land.mask(air)

regional_mean = air.air.groupby(mask).mean()

regional_mean = regional_mean.assign_coords({mask.name: mask.attrs["flag_values"]})
regional_mean = regional_mean.assign_coords({"flags": mask.attrs["flag_meanings"]})

# or
regional_mean = regional_mean.assign_coords({mask.name: mask.attrs["flag_meanings"]})

With our new Grouper API, this should just be a FlagGrouper, it should have all the hooks necessary. We can treat the values as being factorized already (so quite cheap), and then set the output coordinate using flag_meanings. Will need to be careful to not assume sortedness.

How does this look?