Stitching irregular polygons of rasters into a large dataset

I’ve inherited a project which uses raster layers of data over a large area. This data is currently being stored in a large set of tif files, each of which contains the rasters for an irregular polygon, and these polygons tile the spatial extent of our dataset.
I’m trying to find a good method to stitch these together into a single xarray dataset. The resulting object will be much too large to fit into memory, so it will need to be created in an incremental fashion. I’m new to working with these structures, and would greatly appreciate any pointers you folks can provide.

It would be a lot easier if you provided a minimal example, say 4 files that you wanted to stitch