Numerical integration of two-dimensional gradient data is an important step for many slope-measuring optical instruments. However, existing methods are limited by low accuracy or data location restrictions. The zonal integration algorithm in this paper is a generalized process that works with unordered data via Taylor series approximations of finite difference calculations. This method does not require iteration, and all significant steps rely on matrix calculations for a least-squares solution. Simultaneous integration and interpolation is achieved with high accuracy and arbitrary data locations.