We describe a method for reconstructing spatially explicit maps of seasonal palaeoclimate variables from site-based reconstructions. Using a 3D-Variational technique, the method finds the best statistically unbiased, and spatially continuous, estimate of the palaeoclimate anomalies through combining the site-based reconstructions and a prior estimate of the palaeoclimate state. By assuming a set of correlations in the error of the prior, the resulting climate is smoothed both from month to month and from grid cell to grid cell. The amount of smoothing can be controlled through the choice of two length-scale values. The method is applied to a set of reconstructions of the climate of the Last Glacial Maximum (ca. 21,000 years ago, yr BP) for southern Europe derived from pollen data with a prior derived from results from the third phase of the Palaeoclimate Intercomparison Project (PMIP3). We demonstrate how to choose suitable values for the smoothing length scales from the datasets used in the reconstruction.