For studying spectral properties of a nonnormal matrix A is an element of Cnxn, information about its spectrum sigma(A) alone is usually not enough. Effects of perturbations on sigma(A) can be studied by computing epsilon-pseudospectra, i.e. the level sets of the resolvent norm function g(z)=||(zI-A)-1||2. The computation of epsilon-pseudospectra requires determining the smallest singular values sigma min(zI-A) for all z on a portion of the complex plane. In this work, we propose a reduced basis approach to pseudospectra computation, which provides highly accurate estimates of pseudospectra in the region of interest, in particular, for pseudospectra estimates in isolated parts of the spectrum containing few eigenvalues of A. It incorporates the sampled singular vectors of zI - A for different values of z, and implicitly exploits their smoothness properties. It provides rigorous upper and lower bounds for the pseudospectra in the region of interest. In addition, we propose a domain splitting technique for tackling numerically more challenging examples. We present a comparison of our algorithms to several existing approaches on a number of numerical examples, showing that our approach provides significant improvement in terms of computational time.