band_as_variable ignores VRT warping

[alpha-beta-soup]

import rasterio as rio
from rasterio import crs
from rasterio.enums import Resampling
from rasterio.vrt import WarpedVRT
import rioxarray
import xarray as xr

resampling = 'average'
warp_args : dict = {
    'resampling': Resampling._member_map_[resampling],
    'crs': crs.CRS.from_epsg(4326), # Input raster must be converted to WGS84 (4326) for H3 indexing
    'warp_mem_limit': 12000
}
with rio.open(path_to_raster) as src:
    with WarpedVRT(src, src_crs=src.crs, **warp_args) as vrt:
        da : xr.Dataset = rioxarray.open_rasterio(
            vrt,
            band_as_variable=True # <---- !!
        ).chunk(**{'x':'auto','y':'auto'})
        logging.info(da['x'])

array([1831067.831059, 1831077.828418, 1831087.825776, ..., 1856771.040311,
       1856781.03767 , 1856791.035028])

In contast, the same code with band_as_variable=False gives:

array([175.666845, 175.666949, 175.667054, ..., 175.972549, 175.972654,
       175.972758])

The latter has correctly respected the WarpedVRT transformation to EPSG:4326, but the former has retained coordinates in my input projection and seems to ignore the VRT CRS. This is a problem because I need to transform to 4326 prior to performing an H3 index on cell data for later system integration.

Environment Information

python = "^3.10"
gdal = "^3.6.2"
geopandas = "^0.12.2"
h3pandas = "^0.2.3"
rioxarray = "^0.13.4"
dask-geopandas = "^0.3.0"
pyarrow = "^11.0.0"
dask = "^2023.3.0"
click = "^8.1.3"

Linux-5.15.0-41-generic-x86_64-with-glibc2.35

Installed with poetry (pyenv).

[snowman2]

Yes, looks like it happens after the band_as_variable: https://github.com/corteva/rioxarray/blob/33a1be6f1bc3e30cb6d8c1a002d04f439dfb510f/rioxarray/_io.py#L1134

[alpha-beta-soup]

I can work around it by doing a table pivot on the band column; but this does make it more awkward to label bands after their descriptions rather than as integers (which is important for my case).