The purpose of this work is to characterize the performance of different interpolation methods in order to build a global model of ionospheric free electron density, directly derived from the inversions of GNSS occultation data. Indeed, the occultation of GNSS satellites, observed from dual-frequency L-band receivers onboard Low Earth Orbiters (LEOs) is providing valuable ionospheric data since the first mission (GPS/MET) in 1995. This is an excellent scenario to derive, in a simple and precise way, the almost instantaneous distribution of ionospheric electron density in a region of few thousands of kilometers from each occultation event. This can be done, for instance, applying the Improved Abel Transform Inversion technique proposed during the last decade by the UPC authors. In such approach the horizontal gradient of the electron density is taken from the Vertical Total Electron Content gradient as proxy. This assumption increases the accuracy of the electron density retrieval, compared with the Classical Abel Transform Inversion, which assumes spherical symmetry in the occultation region. The availability of numerous occultation events (up to 2000 per day during the full availability of the FORMOSAT-3/COSMIC LEO constellation) enables the execution of studies on the feasibility of simple but accurate ways to interpolate such fields of electron density values, given by each occultation, in an overall electron density model