We consider the classical problem of estimating the covariance matrix of a sub-Gaussian distribution from i.i.d. samples in the novel context of coarse quantization, that is, instead of having full knowledge of the samples, they are quantized to one or two bits per entry. This problem occurs naturally in signal processing applications. We introduce new estimators in two different quantization scenarios and derive nonasymptotic estimation error bounds in terms of the operator norm. In the first scenario, we consider a simple, scale-invariant one-bit quantizer and derive an estimation result for the correlation matrix of a centered Gaussian distribution. In the second scenario, we add random dithering to the quantizer. In this case, we can accurately estimate the full covariance matrix of a general sub-Gaussian distribution by collecting two bits per entry of each sample. In both scenarios, our bounds apply to masked covariance estimation. We demonstrate the near optimality of our error bounds by deriving corresponding (minimax) lower bounds and using numerical simulations.