The theory for the generation of Wannier functions within the generalized Pipek-Mezey approach (Lehtola, S.; Jónsson, H. J. Chem. Theory Comput. 2014, 10, 642) is presented and an implementation thereof is described. Results are shown for systems with periodicity in one, two, and three dimensions as well as isolated molecules. The generalized Pipek-Mezey Wannier functions (PMWF) are highly localized orbitals consistent with chemical intuition where a distinction is maintained between σ- and π-orbitals. The PMWF method is compared with the so-called maximally localized Wannier functions (MLWFs) that are frequently used for the analysis of condensed matter calculations. Whereas PMWFs maximize the localization criterion of Pipek and Mezey, MLWFs maximize that of Foster and Boys and have the disadvantage of mixing σ- and π-orbitals in many cases. The PMWF orbitals turn out to be as localized as the MLWF orbitals as evidenced by cross-comparison of the values of the PMWF and MLWF objective functions for the two types of orbitals. Our implementation in the atomic simulation environment (ASE) is compatible with various representations of the wave function, including real-space grids, plane waves, and linear combinations of atomic orbitals. The projector-augmented wave formalism for the representation of atomic core electrons is also supported. Results of calculations with the GPAW software are described here, but our implementation can also use output from other electronic structure software such as ABINIT, NWChem, and VASP.