The phase of the electronic wave function is not directly measurable but, quite remarkably, it becomes accessible in pairs of isospectral shapes, as recently proposed in the experiment by Moon et al. [Science 319, 782 (2008)]. The method is based on a special property, called transplantation, which relates the eigenfunctions of the isospectral pairs, and allows us to extract the phase distributions, if the amplitude distributions are known. We numerically simulate such a phase extraction procedure in the presence of disorder, which is introduced both as Anderson disorder and as roughness at edges. With disorder, the transplantation can no longer lead to a perfect fit of the wave functions, however we show that a phase can still be extracted-defined as the phase that minimizes the misfit. Interestingly, this extracted phase coincides with (or differs negligibly from) the phase of the disorder-free system, up to a certain disorder amplitude, and a misfit of the wave functions as high as ~5%, proving a robustness of the phase extraction method against disorder. However, if the disorder is increased further, the extracted phase shows a puzzle structure, no longer correlated with the phase of the disorder-free system. A discrete model is used, which is the natural approach for disorder analysis. We provide a proof that discretization preserves isospectrality and the transplantation can be adapted to the discrete systems.