Abstract Stellar atmosphere modelling predicts the luminosity and temperature of a star, together with parameters such as the effective gravity and the metallicity, by reproducing the observed spectral energy distribution. Most observational data come from photometric surveys, using a variety of passbands. We herein present the Python Stellar Spectral Energy Distribution (PySSED) routine, designed to combine photometry from disparate catalogues, fit the luminosity and temperature of stars, and determine departures from stellar atmosphere models such as infrared or ultraviolet excess. We detail the routine’s operation, and present use-cases on both individual stars, stellar populations, and wider regions of the sky. PySSED benefits from fully automated processing, allowing fitting of arbitrarily large data sets at the rate of a few seconds per star.