from osgeo import gdal
import rioxarray
gdal.BuildVRT("merged_mosiac.vrt", ["file1.tiff", "file2.tiff"])
da = rioxarray.open_rasterio("merged_mosiac.vrt")