We study rotating boson star data for numerical relativity using a 3 + 1 decomposition adapted to a curvilinear axisymmetric spacetime with regularization at the rotation axis. The Einstein-Klein-Gordon equations result in a system of six coupled, elliptic, nonlinear equations with an added unknown for the scalar field's frequency omega. Utilizing a Cartesian two-dimensional grid, fourth-order finite differences, and global Newton methods, we generated seven data sets characterized by a rotation azimuthal integer l is an element of [0, 6]. Our numerical implementation is shown to correctly converge both with respect to the resolution size and boundary extension. Thus, global parameters such as the Komar masses and angular momenta can be precisely calculated to characterize these spacetimes. Furthermore, analyzing each family at fixed rotation integer l produces maximum masses and minimum rotation frequencies. Our results coincide with previous results in literature for l is an element of [0, 2] [as in references (Yoshida and Eriguchi 1997 Phys. Rev. D 56 6370; Lai 2004 PhD Thesis; Grandclement, et al 2014 Phys. Rev. D 90 024068; Liebling and Palenzuela 2017 Living Rev. Relativ. 20)], and are new for l > 2. In particular, the study of high-amplitude and localized scalar fields in axial symmetry is revealed to be only possible by adding the sixth regularization variable, thus reaffirming previous work on the importance of regularization in curvilinear-based numerical relativity. These results also provide the groundwork for extending this research to other self-gravitating systems, such as rotating boson stars with nonlinear self-interactions, or the case of massive vector fields.