We consider the classical problem of estimating the covariance matrix of a subgaussian distribution from i.i.d. samples in the novel context of coarse quantization, i.e., 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 non-asymptotic 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 subgaussian 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.