Treating area-level case counts as an aggregated spatial point process offers two advantages over the more common spatially discrete models. First, the spatial correlation can be modelled and interpreted independently of the structure of regions on which data are collected; and second, changes to boundaries over time are easily accommodated. The Expectation-Maximization-Smoothing algorithm has proved adept at estimating intensity surfaces for aggregated point processes, a task for which most other point processes models are unsuitable for. While initially used as a non-parametric method, the EMS algorithm also enables fully model-based inference with a latent-Gaussian intensity surface. This talk will describe this methodology and present newly released software for implementing it.