In ionospheric tomography, the atmospheric electron density is reconstructed from different electron density related measurements, most often from ground-based measurements of satellite signals. Typically, ionospheric tomography suffers from two major complications. First, the information provided by measurements is insufficient and additional information is required to obtain a unique solution. Second, with necessary spatial and temporal resolutions, the problem becomes very high dimensional, and hence, computationally infeasible. With Bayesian framework, the required additional information can be given with prior probability distributions. The approach then provides physically quantifiable probabilistic interpretation for all model variables. Here, Gaussian Markov random fields (GMRFs) are used for constructing the prior electron density distribution. The use of GMRF introduces sparsity to the linear system, making the problem computationally feasible. The method is demonstrated over Fennoscandia with measurements from global navigation satellite system (GNSS) and low Earth orbit (LEO) satellite receiver networks, GNSS occultation receivers, LEO satellite Langmuir probes, and ionosonde and incoherent scatter radar measurements.