Quantitative characterization of discontinuities is fundamental to define the mechanical behavior of discontinuous rock masses. Several techniques for the semi-automatic and automatic extraction of discontinuities and their properties from raw or processed point clouds have been introduced in the literature to overcome the limits of conventional field surveys and improve data accuracy. However, most of these techniques do not allow characterizing flat or subvertical outcrops because planar surfaces are difficult to detect within point clouds in these circumstances, with the drawback of undersampling the data and providing inappropriate results. In this case, 2D analysis on the fracture traces are more appropriate. Nevertheless, to our knowledge, few methods to perform quantitative analyses on discontinuities from orthorectified photos are publicly available and do not provide a complete characterization. We implemented scanline and window sampling methods in a digital environment to characterize rock masses affected by discontinuities perpendicular to the bedding from trace maps, thus exploiting the potentiality of remote sensing techniques for subvertical and low-relief outcrops. The routine, named QDC-2D (Quantitative Discontinuity Characterization, 2D) was compiled in MATLAB by testing a synthetic dataset and a real case study, from which a high-resolution orthophoto was obtained by means of Structure from Motion technique. Starting from a trace map, the routine semi-automatically classifies the discontinuity sets and calculates their mean spacing, frequency, trace length, and persistence. The fracture network is characterized by means of trace length, intensity, and density estimators. The block volume and shape are also estimated by adding information on the third dimension. The results of the 2D analysis agree with the input used to produce the synthetic dataset and with the data collected in the field by means of conventional geostructural and geomechanical techniques, ensuring the procedure’s reliability. The outcomes of the analysis were implemented in a Discrete Fracture Network model to evaluate their applicability for geomechanical modeling.