Traction Force Microscopy (TFM) is a versatile tool to quantify cell-exerted forces by imaging and tracking fiduciary markers embedded in elastic substrates. The computations involved in TFM are ill-conditioned, and data smoothing or regularization is required to avoid overfitting the noise in the tracked substrate displacements. Most TFM calculations depend critically on the heuristic selection of regularization (hyper)parameters affecting the balance between overfitting and smoothing. However, TFM methods rarely estimate or account for measurement errors in substrate deformation to adjust the regularization level accordingly. Moreover, there is a lack of tools to quantify how these errors propagate to the recovered traction stresses. These limitations make it difficult to interpret TFM readouts and hinder comparing different experiments. This manuscript presents an uncertainty-aware TFM technique that estimates the variability in the magnitude and direction of the traction stress vector recovered at each point in space and time of each experiment. In this technique, substrate deformation and its uncertainty are quantified using a non-parametric bootstrap PIV method by resampling the microscopy image pixels (PIV-UQ). This information is passed to a hierarchical Bayesian framework that automatically selects its hyperparameters to perform spatially adaptive regularization conditioned on image quality and propagates the uncertainty to the traction stress readouts (TFM-UQ). We validate the performance of PIV-UQ and TFM-UQ using synthetic datasets with prescribed image quality variations and demonstrate the application of PIV-UQ and TFM-UQ to experimental datasets. These studies show that TFM-UQ locally adapts the level of smoothing, outperforming traditional regularization methods. They also illustrate how uncertainty-aware TFM tools can be used to objectively choose key image analysis parameters like PIV-UQ interrogation window size. We anticipate that these tools will allow for decoupling biological heterogeneity from measurement variability and facilitate automating the analysis of large datasets by parameter-free, input data-based regularization.