We propose an efficient numerical method to obtain local order parameter in two-dimensional systems using spiral boundary conditions. As a benchmark, we first estimate the magnitude of staggered magnetization for the S = 1=2 XXZ Heisenberg model on a square lattice in the whole range of the XXZ anisotropy by density-matrix renormalization group technique. The validity of our method is confirmed by comparing our results with the previous analytical and numerical studies. Then, as further demonstration, we apply our method to obtain the staggered magnetization for the higher-spin cases from S = 1 to 6. The accuracy of the obtained results is validated using the series expansion and the spin-wave theory.