Help me model an implied volatility surface using R or Python

I have quality implied volatility data from 10 days to expiration to 720 days to expiration, across all strikes on a number of underlyings. I would like to understand how to smoothly interpolate these into a neat surface that I can then use to price options using Black-Scholes based on a generic variables (i.e. 50-Delta, 34 days to expiration). I've attached what the raw data looks like hereto.