This research focuses on extracting the statistical features, in space and time, of the monthly rainfall in Saudi Arabia (SA) and the relation to the large-scale atmospheric variability through teleconnection for strategic water resources planning. These features are useful for future predictions. 28 stations distributed over SA for a period between 1970 and 2012 are utilized. According to the Kolmogorov–Smirnov (K-S) test, the Log-normal and Gamma distributions are dominant, while for the Chi-squared (Chi2) test, the Beta distribution is dominant. The K-S is preferable since it works with the original data rather than the Chi² that uses binning, and therefore, some information is lost. The L-moment analysis showed that Person type III is dominant for the wet season while there is no obvious distribution for the dry season. Empirical Orthogonal Function (EOF) analysis is applied to seasonal rainfall data for studying the dominant modes of climate variability and associated large-scale circulation patterns. Our results demonstrate a robust relationship between the wet season (November – April) with El Niño Southern Oscillation (ENSO) and Pacific Decadal Oscillation (PDO), whereas the dry season (June – September) is associated with the Indian Ocean Dipole (IOD). Moreover, the warm (cold) phase of PDO is associated with excess (deficit) rainfall, indicating some predictability of the seasonal rainfall over SA.