Having a map such as:
How to find major directions of variation i.e., for the example given here: north-east--south-west? in an elegant, quick and clean way? Any ideas are more than welcome.
Such a map can be, and usually is, represented in many different ways, including as a grid, a TIN (triangular interpolated network), a set of contours (each of which is given as a contour level together with an ordered sequence of coordinates of points along the contours), a function (specified sometimes as a linear combination of spherical harmonics or other basis functions), an implicit function (such as the zero set of a function), and in other more obscure ways. Usually the map was derived originally from data that can range from spot elevations and topographic features (such as streams and ridgelines) to regular grids of remotely sensed data (as in the Shuttle SRTM elevation datasets) to scanned versions of paper topographic maps (a "DRG"). The commonest representations tend to be grids and sets of contours. (A DRG has to be post-processed into one of the other formats to be useful for this kind of analysis.)
When the map is given as a grid of elevations (a "DEM", or Digital Elevation Model), treat the elevations as if they were unnormalized density values, and compute their first and second moments (as discrete approximations to the integrals: this is a very simple thing to do). From these it is routine to compute the central second moments: the variance-covariance matrix. Its eigenvectors point in the principal directions.
Another approach is to compute the aspect grid for the DEM and post-process that (e.g., by computing average cosine and sine of the aspect). This could be interpreted as an estimate of the minor direction of variation. The aspect is usually computed from a moving-window least-squares fit of a plane or quadratic surface; most formulas relate the aspect to linear functions of a 3 by 3 local neighborhood and therefore can be computed using Fast Fourier Transform techniques, which will be very fast.
When the map is given as a set of contour lines, you can compute a "zonal geometry" for each such closed line. The output will include the major directions of the approximating ellipses. Any grid can produce a set of contour lines, so this approach is available for DEMs, too.
Many GISes provide these calculations as well as methods to convert between various data formats (such as computing contour lines for DEMs).