  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
$S=6$. The accuracy of the obtained results is validated using the series
expansion and the spin-wave theory.
