Pair-instability supernovae (PISNe) have long been predicted to be the final fates of near-zero-metallicity very massive stars ($Z \lt Z_\odot /3$, M$_\mathrm{ZAMS} \gtrsim 140\, \text{M}_\odot$). However, no definite PISN has been observed to date, leaving theoretical modelling validation open. To investigate the observability of these explosive transients, we combine detailed stellar evolution models for PISNe formation, computed from the binary population and spectral synthesis code suite, bpass, with the star formation history of all individual computational elements in the Illustris-TNG simulation. This allows us to compute comic PISN rates and predict their host galaxy properties. Of particular importance is that IllustrisTNG galaxies do not have uniform metallicities throughout, with metal-enriched galaxies often harbouring metal-poor pockets of gas where PISN progenitors may form. Accounting for the chemical inhomogeneities within these galaxies, we find that the peak redshift of PISNe formation is $z=3.5$ instead of the value of $z=6$ when ignoring chemical inhomogeneities within galaxies. Furthermore, the rate increases by an order of magnitude from 1.9 to 29 PISN Gpc$^{-3}$ yr$^{-1}$ at $z=0$, if the chemical inhomogeneities are considered. Using state-of-the-art theoretical PISN light curves, we find an observed rate of 13.8 (1.2) visible PISNe per year for the Euclid-Deep survey, or 83 (7.3) over the 6-yr lifetime of the mission when considering chemically inhomogeneous (homogenous) systems. Interestingly, only 12 per cent of helium PISN progenitors are sufficiently massive to power a superluminous supernova event, which can potentially explain why PISN identification in time-domain surveys remains elusive and progress requires dedicated strategies.