We present a parallel hp-adaptive high order (spectral) discontinuous Galerkin method for approximation of the incompressible Navier-Stokes equations. The spatial discretization consists of equal-order polynomial approximations of the fluid velocity and pressure via discontinuous Galerkin spatial discretizations. For the nonlinear convective term we select the local Lax-Friedrichs flux, while for the divergence and gradient operators central fluxes are chosen. For the diffusive term, we use an interior penalty discontinuous Galerkin method to ensure stability and invertibility. The temporal discretization is an implicit-explicit Runge-Kutta method paired with a high-order splitting procedure to efficiently enforce the incompressibility condition at each time step. The compact stencil size, explicit time stepping of nonlinear terms, and inversion of sparse linear systems make the resulting method simple to parallelize while the local nature of the discontinuous Galerkin approximation makes hp-adaptive refinement natural to implement. We detail our implementation consisting of a tensor product basis of high order polynomials on quadrilateral elements, and implement hp-adaptivity using an inexpensive a posteriori error estimator to determine where refinement is necessary. p-Multigrid and pressure projection techniques are used to precondition the conjugate gradient linear solvers. We present several numerical tests to demonstrate the efficacy of the method, in particular in reducing the number of degrees of freedom needed and allocating computing resources to regions of sharp variation in transient incompressible Navier-Stokes flows.