Many current numerical relativity codes (in three spatial dimensions) share several features: they adopt the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation of Einstein's equations, use finite-difference methods in Cartesian coordinates, and adopt moving puncture coordinates, i.e. ~a combination of 1+log slicing and the "Gamma-driver" condition. While Cartesian coordinates are well suited for many applications, in particular simulations of binaries, spherical polar coordinates have some desirable properties for simulations of single objects, for example gravitational collapse and supernova explosions. We have therefore developed and implemented a new approach that applies in spherical polar coordinates the numerical methods that have previously proven to be extremely successful in Cartesian coordinates. This approach relies on three key ingredients: a reference-metric formulation of the BSSN equations, factoring out appropriate geometrical factors from tensor components, and using a ''partially implicit" Runge-Kutta (PIRK) method. The resulting equations are still singular at the origin of the coordinate system and on the polar axis, but all singular terms can be handled analytically, and the PIRK method is stable even in the presence of these singular terms. Our approach therefore does not rely on a regularization of the equations, and can be used even in the absence of spherical or axi-symmetry. We also applied a reference-metric approach to the formulation of relativistic hydrodynamics, and implemented the resulting equations to perform what we believe are the first self-consistent and stable simulations of general relativistic hydrodynamics in dynamical spacetimes in spherical polar coordinates without the need for regularization or symmetry assumptions.