Modern high-resolution numerical models used in engineering often produce multidimensional maps of outputs (e.g. nodal displacements on a FEM mesh) that may result in more than 10^5 highly correlated outputs for each set of model parameters. Most available meta-modelling techniques, however, are not yet suitable for handling such large maps, including Polynomial Chaos Expansions (PCE). Indeed, the PCE of a numerical model with many outputs is traditionally handled by independently meta-modelling each one of them. We introduce a two-stage PCE approach that aims at solving this problem: in the first stage, PCE is used to compress the map of outputs on a much sparser basis in the map coordinates; in the second stage, standard PCE of the compressed map is carried out w.r.t. the underlying model parameters. Standard PCE post-processing techniques are then used to derive analytical expressions for several stochastic properties of the resulting compressive PCE.