creating geometry variables from bounds

I wonder if it would be in scope for the xvec accessor to have a method that can be used to convert two bounds coordinates to shapely geometry columns?

I'd imagine something like this:

import xarray as xr
import cf_xarray
import xvec

ds = (
    xr.tutorial.open_dataset("air_temperature", chunks={})
    .cf.add_bounds(["lat", "lon"])
    .xvec.bounds_to_shapely(["lat_bounds", "lon_bounds"])
)

That would require xvec to broadcast or stack lat and lon, so we might also require doing that first (but ideally, we'd support 2D geometry variables):

ds = (
    xr.tutorial.open_dataset("air_temperature", chunks={})
    .cf.add_bounds(["lat", "lon"])
    .stack(cells=["lat", "lon"])
    .xvec.bounds_to_shapely(["lat_bounds", "lon_bounds"])
)

Start at the xvec accessor and evaluate the proposed bounds_to_shapely entry point using the examples in the issue. Compare the unstacked and stacked lat/lon cases, including the requirement to broadcast or stack bounds and the possibility of 2D geometry variables. Done means the supported input shape and conversion behavior are specified and implemented.