We perform a new Monte Carlo QCD analysis of pion parton distribution functions, including, for the first time, transverse momentum dependent pion-nucleus Drell-Yan cross sections together with $p_{rm T}$-integrated Drell-Yan and leading neutron electroproduction data from HERA. We assess the sensitivity of the Monte Carlo fits to kinematic cuts, factorization scale, and parametrization choice, and we discuss the impact of the various data sets on the pions quark and gluon distributions. This study provides the necessary step towards the simultaneous analysis of collinear and transverse momentum dependent pion distributions and the determination of the pions 3-dimensional structure.