A method to calculate the crystal field parameters {it ab initio} is proposed and applied to trivalent rare earth impurities in yttrium aluminate and to Tb$^{3+}$ ion in TbAlO$_3$. To determine crystal field parameters local Hamiltonian expressed in basis of Wannier functions is expanded in a series of spherical tensor operators. Wannier functions are obtained by transforming the Bloch functions calculated using the density functional theory based program. The results show that the crystal field is continuously decreasing as the number of $4f$ electrons increases and that the hybridization of $4f$ states with the states of oxygen ligands is important. Theory is confronted with experiment for Nd$^{3+}$ and Er$^{3+}$ ions in YAlO$_3$ and for Tb$^{3+}$ ion in TbAlO$_3$ and a fair agreement is found.