We have inverted velocity solutions from nine geodetic networks distributed across New Zealand to derive present-day continuous horizontal velocity and strain rate fields at the Earth's surface throughout the country. The nine networks contain a total of 362 Global Positioning System (GPS) stations that have been observed at least twice and at least a year apart between 1991 and 1998. The model velocity field is expanded as bicubic spline interpolation functions defined within a curvilinear grid that covers the country and extends into the assumed rigid Australian and Pacific plates to west and east. The inversion jointly minimizes the magnitudes of fitted strain rates and the misfit to the observed velocity data. The spline technique allows high spatial resolution of strain rate variations, especially in regions with spatially dense GPS data. Expansion of the model velocity field on the surface of a sphere allows arbitrarily large areas to be studied. Previously known aspects of New Zealand plate boundary deformation are highlighted with improved resolution, including back-arc extension in the Taupo Volcanic Zone, rotation of the Hikurangi forearc away from this zone in the North Island, and a band of high shear strain rate under the Southern Alps to the southeast of the Alpine fault. Our results reveal several new features, including a region of enhanced shear straining apparently associated with strike-slip faults in the southern North Island and a band of contractional straining subparallel and well east of the Alpine fault that is similar to features found in numerical and sandbox models of continental collision. ¿ 2001 American Geophysical Union