NRPyElliptic sets up initial data (ID) for numerical relativity (NR) using the same numerical methods employed for solving hyperbolic evolution equations. The code implements a hyperbolic relaxation method to solve complex nonlinear elliptic PDEs for NR ID. The hyperbolic PDEs are evolved forward in (pseudo)time, resulting in an exponential relaxation of the arbitrary initial guess to a steady state that coincides with the solution of the elliptic system. The package solves these equations on highly efficient numerical grids exploiting underlying symmetries in the physical scenario. NRPyElliptic is built in the NRPy+ (ascl:1807.025) framework, which facilitates the solution of hyperbolic PDEs on Cartesian-like, spherical-like, cylindrical-like, or bispherical-like numerical grids.