Remote sensing of the earth with satellites yields datasets that can be massive in size, nonstationary in space, and non-Gaussian in distribution. To overcome computational challenges, we use the reduced-rank spatial random effects (SRE) model in a statistical analysis of cloud-mask data from NASA's Moderate Resolution Imaging Spectroradiometer (MODIS) instrument on board NASA's Terra satellite. Parameterisations of cloud processes are the biggest source of uncertainty and sensitivity in different climate models' future projections of Earth's climate. An accurate quantification of the spatial distribution of clouds, as well as a rigorously estimated pixel-scale clear-sky-probability process, is needed to establish reliable estimates of cloud-distributional changes and trends caused by climate change. Here we give a hierarchical spatial-statistical modelling approach for a very large spatial dataset of 2.75 million pixels, corresponding to a granule of MODIS cloud-mask data, and we use spatial change-of-support relationships to estimate cloud fraction at coarser resolutions. Our model is non-Gaussian; it postulates a hidden process for the clear-sky probability that makes use of the SRE model, EM-estimation, and optimal (empirical Bayes) spatial prediction of the clear-sky-probability process. Measures of prediction uncertainty are also given.