High-frequency gravitational waves are the subject of rapidly growing interest in the theoretical and experimental community. In this work we calculate the resonant conversion of gravitational waves into photons in the magnetospheres of neutron stars via the inverse Gertsenshtein mechanism. The resonance occurs in regions where the vacuum birefringence effects cancel the classical plasma contribution to the photon dispersion relation, leading to a massless photon in the medium which becomes kinematically matched to the graviton. We set limits on the amplitude of a possible stochastic background of gravitational waves using X-ray and IR flux measurements of neutron stars. Using Chandra (2–8 keV) and NuSTAR (3–79 keV) observations of RX J1856.6-3754, we set strain limits hclim≃10-26–10-24 in the frequency range 5×1017  Hz≲f≲2×1019  Hz. Our limits are many orders of magnitude stronger than existing constraints from individual neutron stars at the same frequencies. We also use recent JWST observations of the Magnetar 4U 0142+61 in the range 2.7×1013  Hz≲f≲5.9×1013  Hz, setting a limit hclim≃5×10-19. These constraints are in complementary frequency ranges to laboratory searches with CAST, OSQAR and ALPS II. We expect these limits to be improved both in reach and breadth with a more exhaustive use of telescope data across the full spectrum of frequencies and targets.