Gaussian rand om fields (GRF) are a fundamental stochastic model for spatiotemporal data analysis. An essential in gredient of GRF is the covariance function that ch aracterizes the joint Gaussian distribution of the field. Commonly used covariance functions give ri se to fully dense and unstructured covariance matr ices\, for which required calculations are notorio usly expensive to carry out for large data. In thi s work\, we propose a construction of covariance f unctions that result in matrices with a hierarchic al structure. Empowered by matrix algorithms that scale linearly with the matrix dimension\, the hie rarchical structure is proved to be efficient for a variety of random field computations\, including sampling\, kriging\, and likelihood evaluation. S pecifically\, with n scattered sites\, sampling an d likelihood evaluation has an O(n) cost and krigi ng has an O(logn) cost after preprocessing\, parti cularly favorable for the kriging of an extremely large number of sites (e.g.\, predict ing on more sites than observed). We demonstrate comprehensive numerical experiments to show the use of the cons tructed covariance functions and their appealing c omputation time. Numerical examples on a laptop in clude simulated data of size up to one million\, a s well as a climate data product with over two mil lion observations. LOCATION:Seminar Room 1\, Newton Institute CONTACT:info@newton.ac.uk END:VEVENT END:VCALENDAR