The three-dimensional (3D) organization of eukaryotic genomes plays an important role in genome function. While significant progress has been made in deciphering the folding mechanisms of individual chromosomes, the principles of the dynamic large-scale spatial arrangement of all chromosomes inside the nucleus are poorly understood. We use polymer simulations to model the diploid human genome compartmentalization relative to nuclear bodies such as nuclear lamina, nucleoli, and speckles. We show that a self-organization process based on a co-phase separation between chromosomes and nuclear bodies can capture various features of genome organization, including the formation of chromosome territories, phase separation of A/B compartments, and the liquid property of nuclear bodies. The simulated 3D structures quantitatively reproduce both sequencing-based genomic mapping and imaging assays that probe chromatin interaction with nuclear bodies. Importantly, our model captures the heterogeneous distribution of chromosome positioning across cells, while simultaneously producing well-defined distances between active chromatin and nuclear speckles. Such heterogeneity and preciseness of genome organization can coexist due to the non-specificity of phase separation and the slow chromosome dynamics. Together, our work reveals that the co-phase separation provides a robust mechanism for encoding functionally important 3D contacts without requiring thermodynamic equilibration that can be difficult to achieve.