We consider the problem of estimating the intensity map of a spatially random phenomenon over a geographical area observed by a sensor network. The spatial phenomenon of interest is modeled using a Gaussian random field specified by its nonlinear mean and covariance functions. Our proposed algorithm includes two stages: a novel greedy sparse recovery algorithm to estimate the parameters of the mean function, and a spatial interpolation stage using an algorithm called simple kriging. Further, we study the application of the proposed algorithm to radio spectrum cartography, and show that it offers a significant advantage in terms of accuracy of map reconstruction compared to existing methods.