A major objective of modelling geophysical features, biological objects, financial processes and many other irregular surfaces and functions is to develop "shape-preserving" methodologies for smoothly interpolating bivariate data with sudden changes in magnitude or spacing. Shape preservation usually means the elimination of extraneous non-physical oscillations. Classical splines do not preserve shape well in this sense. Empirical experiments have shown that the recently proposed cubic L₁ splines are cable of providing C₁-smooth, shape-preserving, multi-scale interpolation of arbitrary data, including data with abrupt changes in spacing and magnitude, with no need for node adjustment or other user input. However, a theoretic treatment of the bivariate cubic L₁ splines is still lacking. The currently available approximation algorithms are not able to generate the exact coefficients of a bivariate cubic L₁ spline. For theoretical treatment and the algorithm development, we propose to solve bivariate cubic L*#8321; spline problems in a generalized geometric programming framework. Our framework includes a primal problem, a geometric dual problem with a linear objective function and convex cubic constraints, and a linear system for dual-to-primal transformation. We show that bivariate cubic L₁ splines indeed preserve linearity under some mild conditions. Since solving the dual geometric program involves heavy computation, to improve computational efficiency, we further develop three methods for generating bivariate cubic L₁ splines: a tensor-product approach that generates a good approximation for large scale bivariate cubic L₁ splines; a primal-dual interior point method that obtains discretized bivariate cubic L₁ splines robustly for small and medium size problems; and a compressed primal-dual method that efficiently and robustly generates discretized bivariate cubic L₁ splines of large size.