A common problem in the analysis of space-time data is to compress a large dataset in order to extract the underlying trends. Empirical orthogonal function (EOF) analysis is a useful tool for examining both the temporal and the spatial variation in atmospherical and physical processes, and a convenient method of performing this is the Singular Value Decomposition (SVD). Many spatio-temporal models for measurements Z(t,s) at location s at time t, can be written as a sum of a systematic component and a residual component: Z = M + E, where Z, M and E are all T x N matrices. Our approach permits modeling of incomplete data matrices using an EM-like iterative algorithm for the SVD. We model the trend, M, by linear combinations of smooth temporal basis functions derived from left (temporal) singular vectors of Z with the dimension of the model chosen by cross-validation. We further decompose by SVD the spatio-temporal residual matrix E computed as residuals from regressions at each site (column) of the observations on smoothed temporal basis functions. Finally we fit an autoregressive model to the columns (time series) of residuals from the SVD of E. Our aim is to illustrate a simple model to characterize trends and model the variability in large spatio-temporal data matrices. The methodology is demonstrated with 30 years of daily temperature data from Sicily; we obtain a good fit and a compact description of the spatio-temporal variability using just a few smoothed singular vectors.