  Observations of neutron stars, whether in binaries or in isolation, provide
information about the internal structure of the most extreme material objects
in the Universe. In this work, we combine information from recent observations
to place joint constraints on the properties of neutron star matter. We use (i)
lower limits on the maximum mass of neutron stars obtained through radio
observations of heavy pulsars, (ii) constraints on tidal properties inferred
through the gravitational waves neutron star binaries emit as they coalesce,
and (iii) information about neutron stars' masses and radii obtained through
X-ray emission from surface hot spots. In order to combine information from
such distinct messengers while avoiding the kind of modeling systematics
intrinsic to parametric inference schemes, we employ a nonparametric
representation of the neutron-star equation of state based on Gaussian
processes conditioned on nuclear theory models. We find that existing
astronomical observations imply $R_{1.4}=12.32^{+1.09}_{-1.47}\,\mathrm{km}$
for the radius of a $1.4\,M_{\odot}$ neutron star and
$p(2\rho_\mathrm{nuc})=3.8^{+2.7}_{-2.9}\times10^{34}\,\mathrm{dyn}/\mathrm{cm}^2$
for the pressure at twice nuclear saturation density at the 90% credible level.
The upper bounds are driven by the gravitational wave observations, while X-ray
and heavy pulsar observations drive the lower bounds. Additionally, we compute
expected constraints from potential future astronomical observations and find
that they can jointly determine $R_{1.4}$ to ${\cal{O}}(1)\,\mathrm{km}$ and
$p(2\rho_\mathrm{nuc})$ to 80% relative uncertainty in the next five years.
