Apparent motion estimated on satellite data is used for example to compute the wind field in meteorology, and surface currents in oceanography. The satellite images display turbulent fluids with strong rotational patterns at different spatial and temporal scales. This specificity necessitates devising adapted methods, allowing to control the divergence and curl of the retrieved motion field. Vector spline methods are very adapted to that purpose. The vector spline problem is defined as finding a motion field that satisfies a temporal conservation equation at selected control points and that minimizes a regularity constraint in all the image domain. An exact solution of this problem can be found for the 2nd order div-curl regularity constraint. The retrieval of the solution does not require an iterative minimization procedure: a dense matrix must be inverted to compute the spline's coefficients. This matrix unfortunately becomes large and ill-conditioned as the number of control points increases, making the vector spline approach unsuitable for processing large satellite images. This paper presents a method called "Partition of Unity and Optical Flow'' (PUOF), based on a decomposition of the spatial domain: local vector splines are computed in subdomains of the image, then merged using a partition of unity algorithm. The resulting motion field is a good approximation of the exact vector spline solution, and its retrieval is numerically stable and computationally affordable even when processing large data sets, as demonstrated by results obtained on sequences of synthetic and meteorological images.