We investigate the possibility of using proper orthogonal decomposition (POD) in reconstructing complete flow fields from gappy data. The incomplete fields are created from DNS snapshots of flow past a circular cylinder by randomly ommiting data points. We first examine the effectiveness of an existing method and subsequently introduce modifications that make the method robust and lead to the maximum possible resolution at a certain level of spatio-temporal gappiness. We simulate three levels of gappiness at approximately 20%, 50% and 80% in order to investigate the limits of applicability of the new procedure. We find that for the two lower levels of gappiness both the temporal and spatial POD modes can be recovered accurately leading to a very accurate representation of the velocity field. The resulting resolution is improved by more than five times compared to the existing method. However, for 80% gappiness only a few temporal modes are captured accurately while the corresponding spatial modes are noisy. We explain this breakdown of the method in terms of a simple perturbation analysis. This new methodology can be a building block in an effort to develop effective data assimilation techniques in fluid mechanics applications.