We propose a Gaussian process (GP) factor analysis approach for modeling multiple spatio-temporal datasets with non-stationary spatial covariance structure. A novel kernel stick-breaking process based mixture of GPs is proposed to address the problem of non-stationary covariance structure. We also propose a joint GP factor analysis approach for simultaneous modeling of multiple heterogenous spatio-temporal datasets. The performance of the proposed models are demonstrated on the analysis of multi-year unemployment rates of various metropolitan cities in the United States and counties in Michigan. © 2011 IEEE.