Abstract We focus on the discussion of modeling processes that are observed at fixed locations of a region (geostatistics). A standard approach is to assume that the process of interest follows a Gaussian Process with some mean and (valid) covariance functions. It is common to model the covariance function as the product between a variance parameter, and a correlation function which is a function of the Euclidean distance between locations. This implies that the distribution of the process is unchanged when the origin of the index set is translated, and the process is invariant under rotation about the origin; that is the process is stationary and isotropic or homogeneous. However, the assumption of stationarity and isotropy (homogeneity) rarely holds in practice. Commonly, the correlation structures of such processes are influenced by local characteristics resulting in different behaviors in neighborhoods of different spatial locations. We review models that allow for heterogeneous covariance structures and point to some avenues of future research.