We propose a variational order method to explore the ground state phase diagrams of spin-orbit coupled bosons in optical lattices.By properly parameterizing the non-uniform spin order with an ansatz function for each possible phase and comparing their energies which are minimized with respect to the variational parameters,we can effectively identify the lowest energy states and depict the phase diagrams.While the Mott-insulator-superfluid phase boundary is computed by the Gutzwiller method,the spin-ordered phase structures of the deep Mott regime as well as the superfluid regime are studied by the variational order method.Our results agree with those by the classical Monte-Carlo simulations or other numerical methods.For anisotropic spin-orbit coupling,we find some spin-ordered phases in the deep Mott insulator regime that have not been reported before.We also identify the uniform superfluid,the stripe and the checkerboard supersolid phases with exotic spin orders in the superfluid regime.