We present a method for the first principles density functional calculation of relativistic all-electron nuclear magnetic resonance chemical shifts using pseudopotentials. The method is based on the gauge including projector augmented wave approach of Pickard and Mauri [Phys. Rev. B 63, 245101 (2001)]. Relativistic effects are included at the level of the scalar-relativistic zeroth-order regular approximation. The method allows chemical shifts of large, low symmetry structures containing heavy elements to be calculated efficiently. We demonstrate its success for a range of Se and Te containing molecules.