The method of Lin and Huang [Lin and Huang [2004] "Decomposition of incident and reflected higher harmonic waves using four wave gauges," Coast. Eng. 51(5), 395-406.] (LH) is improved for resolution of incident and reflected strongly nonlinear regular waves in shallow waters with measurements of four stationary wave gauges. For first harmonics, wavenumbers, amplitudes and initial phases are obtained by using a nonlinear least squares method. For higher harmonics, wavenumbers of free and bound modes are determined from linear dispersion relation and multiple of first-harmonic wavenumbers, respectively, and the other unknowns are solved by using a linear least squares method. Auto-correlation function is used to determine fundamental wave period for gaining a good performance of Fourier transform. The efficiency and accuracy of the present method are demonstrated by using artificial data and numerical flume data. It is also demonstrated that the present method is less sensitive to signal noise and gauge spacings. Comparison between the present method and the LH method indicates the necessity of employing nonlinear method in determining fundamental wavenumbers of nonlinear shallow-water waves. Finally, the present method is extended to account for obliquely-incident waves. Sensitivity tests indicate the robustness of the extended method with respect to incident angles. Relative position of gauges in the array for avoiding singularity is suggested.