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 105 highly correlated outputs for each set of model parameters. Most available metamodelling 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 metamodelling 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.