Spatial interaction (SI) represents meaningful human relations between areas on the Earth’s surface, such as the reciprocal relations and flows of all kinds among industries, markets, regions, cities, or logistics centers. With the widespread adoption of location-aware technologies and the global diffusion of geographic information systems (GIS), spatial interaction data have been remarkably enriched. In this dissertation research, I develop three unique but closely related exploratory spatial flow data analysis (ESFDA) methods, as an answer to the challenges and opportunities brought by the recent data revolution. Each new method stems from one or more of the following methodological subfields: geovisualization, spatial data mining, and spatial statistics. The first method, dubbed Flow K-function, is a spatial statistical approach to detect spatial clustering patterns of flow data. In other words, it upgrades the classical hot spot detection method to the stage of "hot flow" detection. A set of spatial proximity measures are designed for flow data by integrating endpoint location, length, and direction. The measures can extract both intra-relationships and inter-relationships of flows and serve as the basis of Flow K-function. The second approach, Flow HDBSCAN, is a hierarchical and density-based spatial flow cluster analysis method. Not only can it extract flow clusters from various situations including varying flow densities, lengths, hierarchies, but it also provides an effective way to reveal the potential hierarchical structure of the clusters. The last method is called FlowAMOEBA. It is a data-driven and bottom-up approach for identifying regions of anomalous spatial interactions, based on which it creates a spatial flow weights matrix. It upgrades A Multidirectional Optimum Ecotope-Based Algorithm (AMOEBA) (Aldstadt and Getis 2006) from areal data to spatial flow data through a proper spatial flow neighborhood definition. The method breaks the tradition that spatial interaction data are always collected and modelled between two comparable predefined geographic units, as it delineates the boundaries of anomalous interacting regions regardless of size, shape, scale, or administrative level. The spatial flow weights matrix based on the identified regions can be used to account for network autocorrelation, thus improving confirmatory studies using spatial interaction modeling. These newly developed methods can be utilized individually for data exploration, pattern detection, and hypothesis development. They can also be used jointly to the same application to take advantage of each method. The results of these methods can further be used to form new hypotheses based on explored interesting patterns, to challenge old theories so as to form new ones, to deepen understanding of spatial interaction process, and to improve related confirmatory studies, thus improving related policy-making or problem solving strategies. Three different use cases are presented as to demonstrate the use of each of the methods. The data include a set of motor-vehicle theft and recovery flows, a set of online iPhone transaction flows on the eBay platform, and county-to-county migration flows. Advantages and limitations of each method are tested and discussed thoroughly. Practical usefulness and application implications are also explored and discussed in each scenario.