The structural properties of dense random packings of identical hard spheres (HS) are investigated. The bond order parameter method is used to obtain detailed information on the local structural properties of the system for different packing fractions φ, in the range between φ = 0.53 and φ = 0.72. A new order parameter, based on the cumulative properties of spheres distribution over the rotational invariant w6, is proposed to characterize crystallization of randomly packed HS systems. It is shown that an increase in the packing fraction of the crystallized HS system first results in the transformation of the individual crystalline clusters into the global three-dimensional crystalline structure, which, upon further densification, transforms into alternating planar layers formed by different lattice types.The model of hard spheres (HS) is of fundamental importance in condensed matter physics and material science since it successfully reproduces the essential structural properties of liquids, crystals, glasses, colloidal suspensions and granular media. Packings of HS have also been used in solving important problems of information and optimization theories [1,2]. In this Letter we focus on the structural properties of dense three dimensional (3D) HS systems at different packing fractions φ = π 6 ρσ 3 , where ρ = N/V is the density of N hard spheres in a system volume V and σ is the diameter of the spheres. In particular, we take a large set of packings composed of N = 10 4 identical spheres with periodic boundary conditions, generated using Jodrey-Tory (JT) [5,6] and Lubashevsky-Stillinger (LS) [7] algorithms as described in detail in Refs. [3,4]. The corresponding packing fractions vary in the ranges φ ≃ 0.53 − 0.71 and φ = 0.58 − 0.68 for JT and LS protocols, respectively. Both ranges include the random close packing state (RCP) at φ c ≃ 0.64 (Bernal limit) [8]. The main purpose is to look into the details on how densification of the disordered solid state affects the local and global order of the HS system.To define the local structural properties of the system we use the bond order parameter method [9], which has been widely used in the context of condensed matter physics [9,10] In this method the rotational invariants of rank l of both second q l (i) and third w l (i) order are calculated for each sphere i in the system from the vectors (bonds) connecting its center with the centers of its N nn (i) nearest neighboring spheres: (2) whereY lm (r ij ), Y lm are the spherical harmonics and r ij = r i −r j are vectors connecting centers of spheres i and j. In Eq. (2) l l l m 1 m 2 m 3 are the Wigner 3j-symbols, and the summation in the latter expression is performed over all the indexes m i = −l, ..., l satisfying the condition m 1 + m 2 + m 3 = 0. In 3D the densest possible packings of identical hard spheres are known to be face centered cubic (fcc) and hexagonal close-packed (hcp) lattices (φ ≃ 0.74). To detect these structures we calculate the rotational invariants q 4 , q 6 , and w 6 for each sphere using the ...