This paper shows how modern ideas of scaling can be used to model topography with various morphologies and also to accurately characterize topography over wide ranges of scales. Our argument is divided in two parts. We first survey the main topographic models and show that they are based on convolutions of basic structures (singularities) with noises. Focusing on models with large numbers of degrees of freedom (fractional Brownian motion (fBm), fractional Levy motion (fLm), multifractal fractionally integrated flux (FIF) model), we show that they are distinguished by the type of underlying noise. In addition, realistic models require anisotropic singularities; we show how to generalize the basic isotropic (self-similar) models to anisotropic ones. Using numerical simulations, we display the subtle interplay between statistics, singularity structure and resulting topographic morphology. We show how the existence of anisotropic singularities with highly variable statistics can lead to unwarranted conclusions about scale breaking. We then analyze topographic transects from four Digital Elevation Models (DEMs) which collectively span scales from planetary down to 50 cm (4 orders of magnitude larger than in previous studies) and contain more than 2×108 pixels (a hundred times more data than in previous studies). We use power spectra and multiscaling analysis tools to study the global properties of topography. We show that the isotropic scaling for moments of order =2 holds to within ±45% down to scales \textasciitilde40 m. We also show that the multifractal FIF is easily compatible with the data, while the monofractal fBm and fLm are not. We estimate the universal parameters (a, C1) characterizing the underlying FIF noise to be (1.79, 0.12), where a is the degree of multifractality (0=a=2, 0 means monofractal) and C1 is the degree of sparseness of the surface (0=C1, 0 means space filling). In the same way, we investigate the variation of multifractal parameters between continents, oceans and continental margins. Our analyses show that no significant variation is found for (a, C1) and that the third parameter H, which is a degree of smoothing (higher H means smoother), is variable: our estimates are H=0.46, 0.66, 0.77 for bathymetry, continents and continental margins. An application we developped here is to use (a, C1) values to correct standard spectra of DEMs for multifractal resolution effects.