We present a weak formulation and discretization of the system discovery problem from noisy measurement data. This method of learning differential equations from data replaces point-wise derivative approximations with local integration and improves on the standard SINDy algorithm by orders of magnitude. Linear transformations associated with local integration are used to construct covariance matrices which enforce discovery of parsimonious best-fit models by accurately scaling the error in the residuals during sequentially-thresholded generalized least squares. In the absence of noise, this so-called Weak SINDy framework (WSINDy) is capable of recovering the correct nonlinearities from synthetic data with error in the recovered coefficients falling below the tolerance of the data simulation scheme. As demonstrated by adding white noise directly to the state variables, WSINDy also naturally accounts for measurement noise, with errors in the recovered coefficients scaling proportionally to the signal-to-noise ratio, while significantly reducing the required number of data points and the size of linear systems involved. Altogether, WSINDy combines the ease of implementation of the SINDy algorithm with the natural noise-reduction of integration to arrive at a more robust and user-friendly method of sparse recovery that correctly identifies systems in both small-noise and large-noise regimes. Examples include nonlinear ODEs (Van der Pol Oscillator, Lorenz system) and PDEs (Allen-Cahn, Kuramato-Sivashinsky, Reaction-Diffusion systems) with sharp transitions and/or chaotic behavior.