Construction of “clusters” is in vohelper.py and uses astropy’s SkyCoords and match_catalog_to_sky (asymmetric!).
For three catalogs, we must perform six sky matches to get pairs, then walk the graph to gather the clusters.
This actually is pure astropy and has nothing to do with PyVO as such. As a matter of fact, it is usually smarter to have the remote sides do the cross matches if at all possible.
In this case, since we don’t have a “master catalogue” to match against, that’s actually hard. For smallish crossmatches, the code in vohelper works reasonably well (but it scales horribly when then number of tables increases; use specialised packages when your problem takes that direction).
What’s happening in that code? sky_coords are astropy.SkyCoord instances (in the example code, there’s a function get coordinates_for_table that makes these for essentially arbitrary tables as long as they’re properly marked up).
The code then goes through all pairs of input SkyCoords and uses their match_to_catalog_sky method to generate pairs of indices into these objects that are the closest pairs (that operation isn’t symmetrical, which is why we compute the matches with all permutations).
The remaining code filters out those pairs that are closer than a limit that’s passed in and adds a new pair of rows to be matched to a set. Each row is designated as a pair of table index and row index within that table.
The rest is a graph problem: If you compute the connected subsets of the graph formed in this way, you’ll have all measurements that are crossmatched together and thus, hopefully, correspond to one object.
Sorry for this excursion. Feel free to ignore this.