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 with extreme or unusual observations. A simulation study is proposed in order to compare empirically the calibrated predictive regions with the estimative ones. The proposed method is then applied to the global model assessment of a deterministic model for the prediction of PM10 levels using data from a network of air quality monitoring stations.