We combine both permanent and campaign GPS data from mostly published work to infer a kinematic solution within western North America plate boundary zone. We interpolate 1603 GPS velocity vectors to infer the horizontal velocity gradient tensor field, which is important for both seismic hazards analysis and for constraints on the dynamics. The interpolation algorithm uses continuous bi-cubic splines on the surface of a sphere. The GPS vectors are matched by the model velocity field, in a specified reference frame, in a weighted least- squares inversion. The reference frames of each geodetic study are determined in the inversion procedure. That is, we seek angular velocities, one for each study, that rotate the vectors into one self-consistent frame of reference. The algorithm involves no a priori bias as to the location of fault zones or blocks. The a posteriori variance-covariance matrix of the model strain rates and rotation rates can be evaluated to determine where block behavior is statistically significant, where deformation rates are significant, as well as the significance of the style of inferred strain rates. These model results can be compared with active fault slip observations. The quantitative assessment of the strain rates and rotation rates can be compared with predictions from competing kinematic and dynamic models. Discrimination between such competing models takes into account the a posteriori uncertainties in the kinematic parameters. In cases where discrimination between competing models is not possible, due to large errors in the kinematic solution, then it is possible to investigate where an increased GPS coverage or accuracy is potentially critical in resolving differences between competing models.