Most earthquakes are considered to be caused by stress accumulating in and subsequently releasing from the crust. To extract non-linear and non-stationary earthquake-induced signals associated with stress accumulation, the Hilbert-Huang transform was utilized to filter long-term movements, short-term noise, and frequency-dependent (annual and semi-annual) variations from surface displacements measured by the global positioning system (GPS) in Taiwan. Earthquake-related surface displacements were expressed as horizontal directions (i.e., GPS azimuths) using the north-south and east-west components of residual GPS data to bypass influences resulted from the inhomogeneous nature of the crust. Analytical results showed that the relationships between earthquake occurrence and the aligned GPS azimuth passed the statistical test of the Molchan's error diagram. Aligned GPS azimuths were in agreement with direction of earthquake-related P axes for 81% (26/32) studied events. Areas with the highest paralleling orientations of GPS azimuths appeared around epicenters several days to weeks before earthquake occurrence. Durations from aligned GPS azimuths to earthquake occurrence are roughly proportional to earthquake magnitude. Similar variations of the GPS azimuths were observed in GPS data containing or excluding co-seismic dislocation (i.e., one day before) in the temporal and spatial domain. These suggest that the aligned GPS azimuth could be a promising anomalous phenomenon for studying crustal deformation before earthquakes.