2D Interpolation with periodic boundary conditions
Another function that could work is scipy.ndimage.interpolation.map_coordinates.
It does spline interpolation with periodic boundary conditions.
It does not not directly provide derivatives, but you could calculate them numerically.
I'm running a simulation on a 2D space with periodic boundary conditions.  A continuous function is represented by its values on a grid.  I need to be able to evaluate the function and its gradient at any point in the space.  Fundamentally, this isn't a hard problem -- or to be precise, it's an almost already solved problem.  The function can be interpolated using a cubic spline with scipy.interpolate.RectBivariateSpline.  The reason it's almost solved is that RectBivariateSpline cannot handle periodic boundary conditions, nor can anything else in scipy.interpolate, as far as I can figure out from the documentation.  
Is there a python package that can do this?  If not, can I adapt scipy.interpolate to handle periodic boundary conditions?  For instance, would it be enough to put a border of, say, four grid elements around the entire space and explicitly represent the periodic condition on it?
[ADDENDUM] A little more detail, in case it matters: I am simulating the motion of animals in a chemical gradient.  The continuous function I mentioned above is the concentration of a chemical that they are attracted to.  It changes with time and space according to a straightforward reaction/diffusion equation.  Each animal has an x,y position (which cannot be assumed to be at a grid point). They move up the gradient of attractant.  I'm using periodic boundary conditions as a simple way of imitating an unbounded space.