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.