In this paper, we propose crossing statistics and its generalization, as a new framework to characterize the anisotropy in a 2D field, e.g. height on a surface, extendable to higher dimensions. By measuring $\nu^+$, the number of up-crossing (crossing points with positive slope at a given threshold of height ($\alpha$)), and $N_{tot}$ (the generalized roughness function), it is possible to distinguish the nature of anisotropy, rotational invariance and Gaussianity of any given surface. For the case of anisotropic correlated self- or multi-affine surfaces (even with different correlation lengths in various directions and/or directional scaling exponents), we analytically derive some relations between $\nu^+$ and $N_{tot}$ with corresponding scaling parameters. The method systematically distinguishes the directions of anisotropy, at $3\sigma$ confidence interval using P-value statistics. After applying a typical method in determining the corresponding scaling exponents in identified anisotropic directions, we are able to determine the kind and ratio of correlation length anisotropy. To demonstrate capability and accuracy of the method, as well validity of analytical relations, our proposed measures are calculated on synthetic stochastic rough interfaces and rough interfaces generated from simulation of ion etching. There are good consistencies between analytical and numerical computations. The proposed algorithm can be mounted with a simple software on various instruments for surface analysis and characterization, such as AFM, STM and etc.