How to find the angle between North and horizon for given altaz coords using pyephem?

Question

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?

Answer 1

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

Answer 2

Whether you want to mark the North Celestial Pole with a cross or dot, or put a mark at the north point on the horizon, you are asking a question that is really a question about your camera lens: how does your particular lens, when at the exact focal length that you used when taking the picture, map the curved expanse of sky upon the flat surface of your camera's sensor?

This is a challenge faced not only in astronomy, but by anyone who takes a picture and then wants to use the image later for surveying or spatial computations.

I have heard that professional astronomers used a library called FITS. My impression is that if you explain to the library what kind of lens your camera has and what kind of distortion it produces, that it can tell you the coordinate of each pixel — which should let you find the point of celestial North:

http://fits.gsfc.nasa.gov/fits_libraries.html

Answer 3

just for the first point of getting the right orientation all could be simplified to:

    positionAngle = np.degrees(np.arctan2(CD12, CD11))

as there is the help from NumPy (np) to handle the arctan function in all cases (quadrants) right.

Michel

Edit: if you would like to detect if the image was mirrored: calculate the determinant of the matrix like

mirrored = (CD11 * CD22 - CD12 * CD21) < 0

and check if it's < 0. Than the transformation has mirroring inside.