Jupiter's moon Io is the most volcanically active body in the Solar System, its surface littered with hundreds of active volcanos varying in brightness on different timescales. Studying the spatial distribution and the time variability of these volcanos informs models of the moon's internal structure. Although the surface of Io has been resolved by spacecraft and ground based observatories, there also exists archival data taken by ground based observatories in near-infrared during occultations of Io by Jupiter and other moons. The occultation light curves encode information on the surface features and have a very long time baseline spanning decades which is needed to constrain the long term variability of the volcanos. While these data have mostly been used for targeted studies of particular volcanos, they have never been systematically analyzed as a whole. We use the starry code (https://rodluger.github.io/starry/v1.0.0/) to build a single probabilistic model for all of the observed occultation light curves. It captures the time variable features on Io's surface with accurate uncertainties and allows for inference using Hamiltonian Monte Carlo. I will present preliminary results on our efforts to infer the time dependent map and also explain the data analysis aspect of the project which is highly relevant to the mapping of exoplanet surfaces.