This paper proposes a method to construct well-calibrated frequentist prediction regions, with particular regard to the highest prediction density regions, which may be useful for multivariate spatial prediction. We consider, in particular, Gaussian random fields, and using a calibrating procedure we effectively improve the estimative prediction regions, because the coverage probability turns out to be closer to the target nominal value. Whenever a closed-form expression for the well-calibrated prediction region is not available, we may specify a simple bootstrap-based estimator. Particular attention is dedicated to the associated, improved predictive distribution function, which can be usefully considered for identifying spatial locations ...