Ways of estimating phase shift
The easy to follow but cpu intensive way

If you want something which is obvious how it works and simple to understand, and you don't mind throwing lots of CPU cycles at it (sounds like you don't, if you're using python), then you could just fit sine waves to the two signals and read the phase out from the fit function.

You'd probably want to use scipy.optimise.curve_fit() to fit something like \$V=A \times \sin \left(2 \pi f t+\phi \right)\$ with \$\phi\$ as the only free parameter.  This could converge on several different values of \$\phi\$, so take it modulo \$2\pi\$, then do the same for the current and subtract the two values of \$\phi\$, taking care to add or subtract \$2\pi\$ to get it in the meaningful range.

This will be much, much more CPU intensive than filtering, or doing realtime homodyne detection, but it's easy to follow, only a few lines long, and CPU cycles are cheap. 

If you can't afford those cpu cycles

Then you might want to try generating two quaderature sine waves in software, then using them to do IQ demodulation of the two signals.  Then calculate the phases and subtract.  This will not be as easy to follow, but could be made very fast (even in python) and will be very, very accurate.

Take a chunk of data which is an integer number of cycles (assuming your 50Hz grid is nice and stable, just take 1s of data), and multiply it point-by-point with each of the generated sine waves.  Then integrate each one using numpy.trapz() to get two scalars.  Take these as the arguments to numpy.arctan2(), and voila, you have the phase of that signal relative to the generated sine wave.  Do the same for the the second signal and subtract the two to get the difference. As above, add or subtract \$2\pi\$ to get it in the meaningful range.