We consider a high-dimensional model in which variables are observed over time and on a spatial grid. The model takes the form of a spatio-temporal regression containing time lags and a spatial lag of the dependent variable. Unlike classical spatial autoregressive models, we do not rely on a predetermined spatial interaction matrix but infer all spatial interactions from the data. That is, assuming sparsity, we estimate the spatial and temporal dependence in a fully data-driven way by penalizing a set of Yule-Walker equations. This regularization can be left unstructured but we also propose more customized shrinkage procedures that follow intuitively when observations originate from spatial grids (e.g. satellite images). Finite sample error bounds are derived and estimation consistency is established in an asymptotic framework wherein the sample size and the number of spatial units diverge jointly. A simulation exercise shows strong finite sample performance compared to competing procedures. As an empirical application, we model satellite measured NO2 concentrations in London. Our approach delivers forecast improvements over a competitive benchmark and we discover evidence for strong spatial interactions.