A method for obtaining spatially localized crystalline orbitals starting from delocalized Bloch functions is proposed. The method, that has been implemented in the LCAO CRYSTAL code, is intrinsic and general for nonconducting systems, and provides a set of well localized Wannier functions that can be used for applications that take advantage of their localized character. Examples are given that illustrate the performances and efficiency of the proposed scheme.