python - Implement r.cross / ArcGIS Combine with GDAL

I am using GDAL to process a lot of rasters. I then need to combine all the arrays to create a single output raster that represents all the unique combinations of values. I.E the same thing as r.cross in GRASS or Combine_sa in arcgis.

I'd like to program a similar thing using the GDAL library. What is the algorithm? How would I handle rasters that are too large to fit in RAM? I.e., when chunking them into smaller arrays, how do I ensure that a combination in the nth chunk is given the same value if it appeared in a different chunk?

I'm struggling to come up with anything that will perform well.