© 2016 ACM. With modern focus on remote sensing technology, such as LiDAR, the amount of spatial data, in the form of massive point clouds, has increased dramatically. Furthermore, repeated surveys of the same areas are becoming more common. This trend will only increase as topographic changes prompt surveys over already scanned areas, in which case we obtain large spatiotemporal datasets. An initial step in the analysis of such spatial data is to create a digital elevation model representing the terrain, possibly over time. In the case of spatial (spatiotemporal, respectively) datasets, these models often represent elevation on a 2D (3D, respectively) grid. This involves interpolating the elevation of LiDAR points on these grid points. In this article, we show how to efficiently perform natural neighbor interpolation over a 2D and 3D grid. Using a graphics processing unit (GPU), we describe different algorithms to attain speed and GPU-memory tradeoffs. Our experimental results demonstrate that our algorithms not only are significantly faster than earlier ones but also scale to much bigger datasets that previous algorithms were unable to handle.