We developed a multivariate Gaussian process model to jointly analyze high-resolution remote sensing data and climate model output. With a basis function representation, the resulting model can achieve efficient computation and to describe potentially non-stationary spatial dependence. The predictive distribution provides statistical downscaling from the coarse-resolution climate model output, borrowing strength spatially and across high-resolution remote sensing data. We implement the proposed method for downscaling Sea Surface Temperature (SST) over the Great Barrier Reef (GBR). Our method reduces the mean squared predictive error by about 20% compared with the state of the art and produces a predictive distribution enabling holistic uncertainty quantification analyses.