Efficient ab initio computational methods for the calculation of thermoelectric transport properties of materials are of great avail for energy harvesting technologies. The BoltzTraP code has been largely used to efficiently calculate thermoelectric coefficients. However, its current version that is publicly available is based only on the constant relaxation time (RT) approximation, which usually does not hold for real materials. Here, we extended the implementation of the BoltzTraP code by incorporating realistic k-dependent RT models of the temperature dependence of the main scattering processes, namely, screened polar and nonpolar scattering by optical phonons, scattering by acoustic phonons, and scattering by ionized impurities with screening. Our RT models are based on a smooth Fourier interpolation of Kohn-Sham eigenvalues and its derivatives, taking into account non-parabolicity (beyond the parabolic or Kane models), degeneracy and multiplicity of the energy bands on the same footing, within very low computational cost. In order to test our methodology, we calculated the anisotropic thermoelectric transport properties of low temperature phase (Pnma) of intrinsic p-type and hole-doped tin selenide (SnSe). Our results are in quantitative agreement with experimental data, regarding the evolution of the anisotropic thermoelectric coefficients with both temperature and chemical potential. Hence, from this picture, we also obtained the evolution and understanding of the main scattering processes of the overall thermoelectric transport in p-type SnSe.