Constraining numerical geodynamo models with surface geomagnetic observations is very important in many respects: it directly helps to improve numerical geodynamo models, and expands their geophysical applications beyond geomagnetism. A successful approach to integrate observations with numerical models is data assimilation, in which Bayesian algorithms are used to combine observational data with model outputs, so that the modified solutions can then be used as initial conditions for forecasts of future physical states. In this paper, we present the first geomagnetic data assimilation framework, which comprises the MoSST core dynamics model, a newly developed data assimilation component (based on ensemble covariance estimation and optimal interpolation), and geomagnetic field models based on paleo, archeo, historical and modern geomagnetic data. The overall architecture, mathematical formulation, numerical algorithms and computational techniques of the framework are discussed. Initial results with 100-year geomagnetic data assimilation and with synthetic data assimilation are presented to demonstrate the operation of the system.