How to calculate terrestial coordinates from celestial coordinates?

[mdhoney]

I have a picture of the night sky and the matching RA/Dec coordinates of the center. I also have the azimuth and elevation the camera was pointed at, and the time. How can I use this information to calculate the latitude and longitude on Earth where the picture was taken from? I've tried working through some of the examples on the Coordinate Transformations page but I can't figure out how to apply that to my problem. Would appreciate any pointers.

[drbitboy]

Long version

Start at the NAIF/SPICE required reading for Rotations:

• Rotations
• Notation
• References

It will take a while, but once you understand it

• the world will never be the same.
• you will not need the rest of the explanation below, because you can figure it out on your own

Short version

There are four 3-D frames (orientations): image (camera); local (latitude and longitude); Earth Body-Fixed (rotation at time); target (RA, DEC). You have the rotations (connection) between pairs of three of the four frames, and can calculate the fourth, local, by difference.

Convert everything to rotation matrices, which matrices convert vectors from one frame to another frame:

• RA, DEC, Twist => Mj2k2img = [J2000]-to-[Image] matrix
• Azimuth, Elevation => Mlocal2img = [Latitude-Longitude]-to-[Image]
• Time => Mj2k2ebf = [J2000]-to-[Earth-Body-Fixed]

We will solve for the unknown rotation, Mebf2local = [Earth-Body-Fixed]-to-[Latitude-Longitude]

Then, since

Mj2k2img = Mlocal2img * Mebf2local * Mj2k2ebf

we can solve for Mebf2local like this

Mebf2local = Mlocal2imgT * Mj2k2img * Mj2k2ebfT

where the [T] suffix means transpose. With the Python interface to the NAIF/SPICE toolkit, that solution it's one line of code.

With another few lines to get from the known angles to the known matrices, and then one or two to get from the solved-for matrix to latitude and longitude, there will be more lines parsing the inputs and writing the output than you will need to solve the problem.

I think you can get away without knowing Twist (Direction of J2000 North, clockwise from up, in the image).

Then convert the solved-for Mebf2local back to the local angles of latitude and longitude.

is azimuth something like "degrees clockwise from north?"
is elevation something like "degrees above horizon-at-azimuth?"
What is the twist e.g. twist is 0 if the plane containing local "up" (-gravity) and "up" in the picture also contains the boresight pointing at RA/DEC?

[mdhoney]

Thanks for the tip. Yes, azimuth is measured clockwise from North (0 degrees) and elevation is measured as degrees above the horizon.
I was hoping to be able to solve this problem within PyEphem though, without having to set up a bunch of rotation matrices. Do you think that will be possible, @drbitboy?

[drbitboy]

Not directly in Pyephem, no, I don't think so. That page does say, if you have RA/DEC and Lat/Lon, that you can calculate alt/az. So you could probably run an algorithm (binary search?) to solve for Lat/Lon, although it might get messy near the PM.

Also, since the RA/DEC frame is essentially co-polar with the Lat/Lon frame at any epoch, Pyephem might allow you to calculate at RA/DEC from the alt/az at a random Lat/Lon (e.g. 0/0) of the observation epoch, then you could adjust the Lat by the difference between that calculated DEC and the known observation DEC, then calculate a second RA/DEC from the alt/az at that adjusted Lat and the same Lon (0), then adjust the Lon by the difference between the second calculated RA and the know observation RA. For higher accuracy that process could probably be repeated by successive substitution or even Newton's method.