I played around with the position_at() function, and managed to get it to work, so that it gave me my landing predictions as a position vector. But... it wasn't accurate enough for me. So I started from scratch on a predictor that took into account drag. Using some equations, I was able to start with a velocity, and simulate what would happen using a basic physics model until it hit the ground. This turned out to be a lot simpler than I thought it would be - the hardest part was figuring out the atmospheric density on Kerbin at different altitudes (because its modeled on the US standard atmosphere).
I wrote some functions (https://github.com/MusicalHQ/KSP-Landing-Predictions/blob/master/Landing.py) for this.
So far, you can use the functions to calculate landing coordinates (right now just position vectors for body.reference_frame) on bodies with no atmosphere, and on Kerbin. I'll try and add Eve, Duna, Laythe, Earth, Venus, Mars, Titan and maybe some OPM moons too.
My trajectories line up pretty well to mechjeb's, but that's the problem. My only way of checking whether my math is correct is comparing my results to mechjeb. If someone could have a look at the code and tell me whether my math makes sense, it would be appreciated.
The only thing holding back the speed of the calculations is drawing the trajectory. I have a few ideas about how to fix this, but right now if you want to calculate the landing site, turn off drawing.