I've got a set of many astrophotos taken with the camera on a tripod. Used a bubble level to make sure that the long side of the frame is parallel to the horizon, and I know the alt/az (and equatorial) coordinates of the center of each photo.
Now I'm writing some python code to overlay an indicator over each image to mark the North direction. Can I use pyephem to get the angle between the North Celestial Pole direction and the horizontal direction for given alt/az coordinates? Any clue?
I would try to create a "World Coordinate System" (WCS) for each of your images, this is essentially a mapping between your pixel coordinates and sky coordinates (i.e. RA and Dec). You can use a tools such as one available at http://astrometry.net to automatically solve your images based on the star patterns visible in the image. This will generate a WCS for the image.
The astrometry.net solver (if you have the appropriate dependancies installed) can generate a png version of your image with known celestial objects marked. This might be enough for your purposes, but if not, you could read the image WCS in to python using the astropy.wcs package and use that to determine the orientation of the image and then mark up the image however you like.
Here's some quick and dirty code which you might try to adapt to your purpose:
import math
import subprocess
import astropy.units as u
import astropy.fits as fits
## Solve the image using the astrometry.net solve-field tool.
## You'll want to look over the options for solve-field and adapt this call
## to your images.
output = subprocess.check_output(['solve-field', filename])
## Read Header of Image (assumes you are working off a fits file with a WCS)
## If not, you can probably read the text header output my astrometry.net in
## a similar fashion.
hdulist = fits.open(solvedFilename)
header = hdulist[0].header
hdulist.close()
CD11 = float(header['CD1_1'])
CD12 = float(header['CD1_2'])
CD21 = float(header['CD2_1'])
CD22 = float(header['CD2_2'])
## This is my code to interpet the CD matrix in the WCS and determine the
## image orientation (position angle) and flip status. I've used it and it
## seems to work, but there are some edge cases which are untested, so it
## might fail in those cases.
## Note: I'm using astropy units below, you can strip those out if you keep
## track of degrees and radians manually.
if (abs(CD21) > abs(CD22)) and (CD21 >= 0):
North = "Right"
positionAngle = 270.*u.deg + math.degrees(math.atan(CD22/CD21))*u.deg
elif (abs(CD21) > abs(CD22)) and (CD21 < 0):
North = "Left"
positionAngle = 90.*u.deg + math.degrees(math.atan(CD22/CD21))*u.deg
elif (abs(CD21) < abs(CD22)) and (CD22 >= 0):
North = "Up"
positionAngle = 0.*u.deg + math.degrees(math.atan(CD21/CD22))*u.deg
elif (abs(CD21) < abs(CD22)) and (CD22 < 0):
North = "Down"
positionAngle = 180.*u.deg + math.degrees(math.atan(CD21/CD22))*u.deg
if (abs(CD11) > abs(CD12)) and (CD11 > 0): East = "Right"
if (abs(CD11) > abs(CD12)) and (CD11 < 0): East = "Left"
if (abs(CD11) < abs(CD12)) and (CD12 > 0): East = "Up"
if (abs(CD11) < abs(CD12)) and (CD12 < 0): East = "Down"
if North == "Up" and East == "Left": imageFlipped = False
if North == "Up" and East == "Right": imageFlipped = True
if North == "Down" and East == "Left": imageFlipped = True
if North == "Down" and East == "Right": imageFlipped = False
if North == "Right" and East == "Up": imageFlipped = False
if North == "Right" and East == "Down": imageFlipped = True
if North == "Left" and East == "Up": imageFlipped = True
if North == "Left" and East == "Down": imageFlipped = False
print("Position angle of WCS is {0:.1f} degrees.".format(positionAngle.to(u.deg).value))
print("Image orientation is North {0}, East {1}.".format(North, East))
if imageFlipped:
print("Image is mirrored.")
## Now you have position angle and flip status and can mark up your image