The current monitoring of modern distribution networks by operators with the aim to improve the network visibility and thus detect and mitigate challenges such as power quality issues has resulted in massive data generation. Such data need to be analysed and technical observations realised. Thus, this paper presents bad data detection for spatiotemporal harmonic power quality field measurements obtained from six power quality analysers (installed at the connection points of customer-owned photovoltaic and small-scale energy storage and the mains units) with high level of granularity (1s interval). A total of 4.1 billion measured values have been evaluated in the bad data detection. Multivariate statistical methods (principal component analysis and Hotelling’s T2 statistic) have been combined with Bag of Little Bootstrap (BLB) to analyse the massive data. Specifically, the BLB technique has been deployed for calculating the upper control limit for the bad data detection, thus, the spatiotemporal nature of the data has been fully captured. The method has been applied to harmonic load flow analysis of a practical distribution network. The findings reveal how the neglect of spatiotemporality could conceal the actual harmonic penetration levels in a data-driven harmonic load flow analysis. Other interesting technical findings are presented.