Merge pull request #110 from CTetford/develop

AutoPA

Author: ClutchplateDude

Software/Addons/AutoPA/autopa_modules/indi.py

@@ -6,7 +6,7 @@
 import io
 from astropy.coordinates import Angle
 from astropy import units as u
-
+    
 class IndiClient(PyIndi.BaseClient):
     def __init__(self):
         super(IndiClient, self).__init__()
@@ -34,7 +34,7 @@ def serverConnected(self):
         pass
     def serverDisconnected(self, code):
         pass
-
+    
 def capture(indiclient, ccd, exposure, outputfile, blobEvent):
     # Let's take some pictures
     device_ccd=indiclient.getDevice(ccd)

Software/Addons/AutoPA/autopa_modules/platesolve.py

@@ -1,11 +1,12 @@
 #!/usr/bin/env python3
 
 import subprocess
-import re, math, os
+import re, math, os, time
 from pathlib import Path
 from astropy.time import Time
 from astropy.coordinates import SkyCoord, EarthLocation, AltAz, Angle
 from astropy import units as u
+from astropy.wcs import WCS
 from autopa_modules import LX200
 import contextlib
 
@@ -38,17 +39,11 @@ def solve_ASTAP (filename, RA, DEC, radius, pixel_resolution):
     cmd = f"/opt/astap/astap -r 30 -speed auto -o /tmp/solution -f {filename}"    
     try:
         output = str(subprocess.check_output(cmd, shell=True))
+        output = WCS("/tmp/solution.wcs")
+        object = SkyCoord(output.wcs.crval[0], output.wcs.crval[1], frame='fk5', unit='deg')
     except:
-        return ("Failed", "Failed")
-    textfile = open("/tmp/solution.wcs", 'r')
-    filetext = textfile.read()
-    textfile.close()
-    try:
-        RA = re.search("CRVAL1\s*=\s*(\d\.\d+E\+\d+)", filetext).group(1)
-        DEC = re.search("CRVAL2\s*=\s*(\d\.\d+E\+\d+)", filetext).group(1)
-    except:
-        return ("Failed", "Failed")
-    return (Angle(RA, unit=u.deg).to_string(unit=u.hour), Angle(DEC, unit=u.deg).to_string(unit=u.degree))
+        return ("Failed")
+    return (object)
 
 def solve_ASTROMETRY (filename, RA, DEC, radius, pixel_resolution):
     config_file="/home/astroberry/.local/share/kstars/astrometry/astrometry.cfg"
@@ -71,7 +66,7 @@ def solve_ASTROMETRY (filename, RA, DEC, radius, pixel_resolution):
         DEC_S = result.group("DEC_S")
         return (f"{RA_H}h{RA_M}m{RA_S}s", f"{DEC_D}d{DEC_M}m{DEC_S}s")
 
-def polarCalc(mylat, mylong, myelev, observing_time, p1RA, p1DEC, p2RA, p2DEC, p3RA, p3DEC):
+def polarCalc(mylat, mylong, myelev, observing_time, p1, p2, p3):
     #iers.conf.auto_download = False
     #iers.conf.auto_max_age = None
 
@@ -82,17 +77,17 @@ def polarCalc(mylat, mylong, myelev, observing_time, p1RA, p1DEC, p2RA, p2DEC, p
     observing_location = EarthLocation(lat=mylat*u.deg, lon=mylong*u.deg, height=myelev*u.m)
 
     #Create coordinate objects for each point
-    p1 = SkyCoord(p1RA, p1DEC, unit='deg')
-    p2 = SkyCoord(p2RA, p2DEC, unit='deg')
-    p3 = SkyCoord(p3RA, p3DEC, unit='deg')
+    #p1 = SkyCoord(p1RA, p1DEC, unit='deg')
+    #p2 = SkyCoord(p2RA, p2DEC, unit='deg')
+    #p3 = SkyCoord(p3RA, p3DEC, unit='deg')
     p1X = (90 - p1.dec.degree) * math.cos(p1.ra.radian)
     p1Y = (90 - p1.dec.degree) * math.sin(p1.ra.radian)
     p2X = (90 - p2.dec.degree) * math.cos(p2.ra.radian)
     p2Y = (90 - p2.dec.degree) * math.sin(p2.ra.radian)
     p3X = (90 - p3.dec.degree) * math.cos(p3.ra.radian)
     p3Y = (90 - p3.dec.degree) * math.sin(p3.ra.radian)
 
-    #Calculate center of circle using three points in the complex plane. DEC is treated as unitless for the purposes of the calculation.
+    #Calculate center of circle using three points in the complex plane. RA/DEC are treated as unitless for the purposes of the calculation.
     x, y, z = complex(p1X,p1Y), complex(p2X,p2Y), complex(p3X,p3Y)
     w = z-x
     w /= y-x
@@ -104,26 +99,30 @@ def polarCalc(mylat, mylong, myelev, observing_time, p1RA, p1DEC, p2RA, p2DEC, p
     resultDEC = (90 - math.sqrt(resultX**2 + resultY**2))
     resultRA = math.atan2(resultY, resultX)*360 / (2*math.pi)
     if resultRA < 0:
-            resultRA = (180-abs(resultRA))+180
+        resultRA = (180-abs(resultRA))+180
 
-    #Create coordinate object for current alignment offset
-    offset = SkyCoord(resultRA, resultDEC, frame='itrs', unit='deg', representation_type='spherical', obstime=observing_time)
     print(f"Current alignment in RA/DEC: {Angle(resultRA*u.deg).to_string(u.hour, precision=2)}/{Angle(resultDEC*u.deg).to_string(u.degree, precision=2)}.")
-
-    #Create coordinate object for pole
-    pole = SkyCoord(0, 90, frame='itrs', unit='deg', representation_type='spherical', obstime=observing_time)
 
-    #Create coordinate object for pole
-    poleAzAlt = pole.transform_to(AltAz(obstime=observing_time,location=observing_location))
-    print(f"True polar alignment in Az./Alt.: 0h00m00s/{poleAzAlt.alt.to_string(u.degree, precision=2)}.")
-
-    #Transform current alignment to Alt/Az coordinate system
-    offsetAzAlt = offset.transform_to(AltAz(obstime=observing_time,location=observing_location))
+    #Create alt/az coordinate object for pole
+    poleAzAlt = RADECtoAltAz(observing_location, observing_time, 0, 90)
+    print(f"True polar alignment in Az./Alt.: 00h00m00s/{poleAzAlt.alt.to_string(u.degree, precision=2)}.")
+    
+    #Create alt/az coordinate object for current alignment
+    offsetAzAlt = RADECtoAltAz(observing_location, observing_time, resultRA, resultDEC)
     print(f"Current alignment in Az./Alt.: {offsetAzAlt.az.to_string(u.hour, precision=2)}/{offsetAzAlt.alt.to_string(u.degree, precision=2)}.")
 
     #Calculate offset deltas from pole
     #Normalize the azimuth values to between -180 and 180 degrees prior to determining offset.
     errorAz = (((poleAzAlt.az.deg + 180) % 360 - 180)-((offsetAzAlt.az.deg + 180) % 360 - 180))*60
     errorAlt = (poleAzAlt.alt.deg-offsetAzAlt.alt.deg)*60
 
-    return errorAz, errorAlt
+    return errorAz, errorAlt
+    
+def RADECtoAltAz(observing_location, observing_time, RA, DEC):
+    #Create coordinate object for current alignment offset
+    offset = SkyCoord(RA, DEC, frame='itrs', unit='deg', representation_type='spherical', obstime=observing_time)
+    
+    #Transform current alignment to Alt/Az coordinate system
+    AzAlt = offset.transform_to(AltAz(obstime=observing_time,location=observing_location))
+
+    return AzAlt

Software/Addons/AutoPA/polaralign_autocalibration.py

@@ -22,7 +22,7 @@ def polarcalibrate(axis, indiclient, iteration, observing_location, telescope, c
         calibration_distance = calibration_distance * -1
 
     #Initialize steppers to a known backlash position
-    initSteppers(indiclient, telescope, serialport, backlashCorrection)
+    initSteppers(axis, indiclient, telescope, serialport, backlashCorrection)
 
     p1RA, p1DEC, p1Time = captureSolve(indiclient, RA, DEC, radius, exposure, telescope, ccd, blobEvent, pixel_resolution, "altazCalibration1.fits")
     p1Time = datetime.utcfromtimestamp(p1Time)
@@ -85,21 +85,21 @@ def captureSolve(indiclient, RA, DEC, radius, exposure, telescope, ccd, blobEven
                 print ("Image could not be solved after 5 attempts. Exiting program.")
                 exit()
         solveTimeFinish = time.time()
-    print (f"Image solved in {(solveTimeFinish - solveTimeStart):.2f} seconds. Coordinates are RA {resultRA}, DEC {resultDEC}")
+    print (f"Image solved in {(solveTimeFinish - solveTimeStart):.2f} seconds. Coordinates are RA {Angle(resultRA, unit=u.deg).to_string(unit=u.hour)}, DEC {Angle(resultDEC, unit=u.deg).to_string(unit=u.degree)}")
 
     return (resultRA, resultDEC, captureTime)
 
-def initSteppers(indiclient, telescope, serialport, backlashCorrection):
+def initSteppers(axis, indiclient, telescope, serialport, backlashCorrection):
 
     print ("Initializing steppers")
     #Disconnect telescope mount from INDI to free up serial port for Alt/Az adjustments
     indi.disconnectScope(indiclient, telescope)
 
     #Increase both Az & Alt by approximately 16-steps
-    LX200.sendCommand(f":MAZ{backlashCorrection[0]}#", serialport)
-    time.sleep(0.5)
-    LX200.sendCommand(f":MAL{backlashCorrection[1]}#", serialport)
-
+    if axis == "az":
+        LX200.sendCommand(f":MAZ{backlashCorrection[0]}#", serialport)
+    elif axis == "alt":
+        LX200.sendCommand(f":MAL{backlashCorrection[1]}#", serialport)
     time.sleep(3)
 
     #Re-connect telescope mount to INDI before disconnecting from the INDI server
@@ -188,11 +188,11 @@ def waitForMotors(serialport):
 else:
     backlashCorrection = [-4.2665,-0.1647]
 
-iteration = input("How many calibration iterations to perform? (Default is 3) ")
+iteration = input("How many calibration iterations to perform? (Default is 4) ")
 try:
     iteration = int(iteration)
 except:
-    iteration = 3
+    iteration = 4
 
 #Create location object based on lat/long/elev
 observing_location = EarthLocation(lat=mylat*u.deg, lon=mylong*u.deg, height=myelev*u.m)

Software/Addons/AutoPA/polaralign_automatic.py

@@ -0,0 +1,169 @@
+#!/usr/bin/env python3
+from astropy import units as u
+from astropy.coordinates import SkyCoord, EarthLocation, AltAz, Angle
+from astropy.time import Time
+from astropy.utils import iers
+import sys, subprocess, time, math
+from autopa_modules import indi
+from autopa_modules import platesolve
+from autopa_modules import LX200
+from datetime import datetime
+from statistics import mean
+import argparse
+
+def slewCaptureSolve(indiclient, RA, DEC, radius, exposure, telescope, ccd, blobEvent, pixel_resolution, filename="platesolve.fits"):
+    print (f"Slewing \"{telescope}\" to {Angle(RA*u.deg).to_string(unit=u.hour)}, {Angle(DEC*u.deg).to_string(unit=u.degree)}")
+    indi.slewSync(indiclient, telescope, RA, DEC)
+    #Wait one second to let the telescope stop moving and let any motor vibrations reduce
+    time.sleep(1)
+    
+    original_exposure = exposure
+    solved = False
+    while not solved:
+        print (f"Capturing {exposure} second exposure on \"{ccd}\"")
+        captureTime = indi.capture(indiclient, ccd, exposure, filename, blobEvent)
+        
+        print (f"Attempting plate solving in the region of: RA {Angle(RA*u.deg).to_string(unit=u.hour)}, DEC {Angle(DEC*u.deg).to_string(unit=u.degree)}.")
+        solveTimeStart = time.time()
+        #resultRA, resultDEC = platesolve.solve_ASTROMETRY(filename, RA, DEC, radius, pixel_resolution)
+        result = platesolve.solve_ASTAP(filename, RA, DEC, radius, pixel_resolution)
+        if not isinstance(result, str):
+            solved = True
+        else:
+            print("Image could not be solved. Attempting another capture with increased exposure time (Max is 10 second exposure)")
+            if exposure < 10:
+                exposure += 1
+            elif exposure - original_exposure >= 5:
+                print ("Image could not be solved. Exiting.")
+        solveTimeFinish = time.time()
+    print (f"Image solved in {(solveTimeFinish - solveTimeStart):.2f} seconds.")
+    #print (f"Image solved. Coordinates are RA {Angle(resultRA, unit=u.deg).to_string(unit=u.hour)}, DEC {Angle(resultDEC, unit=u.deg).to_string(unit=u.degree)}")
+
+    return (result, captureTime)
+    
+def initSteppers(indiclient, telescope, serialport, backlashCorrection):
+
+    print ("Initializing steppers")
+    #Disconnect telescope mount from INDI to free up serial port for Alt/Az adjustments
+    indi.disconnectScope(indiclient, telescope)
+    
+    #Increase both Az & Alt by approximately 16-steps
+    LX200.sendCommand(f":MAZ{backlashCorrection[0]}#", serialport)
+    LX200.sendCommand(f":MAL{backlashCorrection[1]}#", serialport)
+
+    time.sleep(3)
+    
+    #Re-connect telescope mount to INDI before disconnecting from the INDI server
+    indi.connectScope(indiclient, telescope)
+    return
+
+parser = argparse.ArgumentParser(usage='%(prog)s [mylat] [mylong] [myelev] [options]', description='OpenAstroTracker AutoPA: This tool is used to automatically rotate the mount, capture images,\
+        plate solve, and calculate the polar alignment error of the OAT. Serial commands are automatically issued to the OAT to adjust the motorized altitude/azimuth axis to correct this error.')
+parser.add_argument("mylat", help="your latitude in degrees", type=float)
+parser.add_argument("mylong", help="your longitude in degrees", type=float)
+parser.add_argument("myelev", help="your elevation in metres", type=float)
+parser.add_argument("--serialport", help="serial port address for the OAT (default is /dev/ttyACM0)", type=str)
+parser.add_argument("--targetRA", help="initial starting RA in degrees (default is 0)", type=float)
+parser.add_argument("--targetDEC", help="initial starting DEC in degrees (default is 85)", type=float)
+parser.add_argument("--exposure", help="exposure time in seconds (default is 8 seconds)", type=float)
+parser.add_argument("--telescope", help="name of INDI telescope to control movement (default is Ekos \"Telescope Simulator\")", type=str)
+parser.add_argument("--ccd", help="name of INDI CCD for capturing exposures (default is Ekos \"CCD Simulator\")", type=str)
+parser.add_argument("--radius", help="field radius of plate solving", type=float)
+parser.add_argument("--pixelSize", help="CCD pixel size in micrometres. Used to decrease plate solving time.", type=float)
+parser.add_argument("--pixelX", help="Quantity of pixels in the X direction.", type=float)
+parser.add_argument("--pixelY", help="Quantity of pixels in the Y direction.", type=float)
+parser.add_argument("--focalLength", help="Lens focal length in millimetres. Used to decrease plate solving time.", type=float)
+parser.add_argument("--nomove", help="Run AutoPA sequence but do not move the steppers.", action="store_true")
+args = parser.parse_args()
+mylat = args.mylat
+mylong = args.mylong
+myelev = args.myelev
+if args.serialport:
+    serialport = args.serialport
+else:
+	serialport = "/dev/ttyACM0"
+if args.targetRA:
+    targetRA = args.targetRA
+else:
+	targetRA = 0
+if args.targetDEC:
+    targetDEC = args.targetDEC
+else:
+	targetDEC = 85
+if args.exposure:
+    exposure = args.exposure
+else:
+	exposure = 8.0
+if args.telescope:
+    telescope = args.telescope
+else:
+	telescope = "Telescope Simulator"
+if args.ccd:
+    ccd = args.ccd
+else:
+	ccd = "CCD Simulator"
+if args.radius:
+    radius = args.radius
+else:
+	radius = 30
+if args.pixelSize and args.focalLength:
+    pixel_resolution = (args.pixelSize/args.focalLength)*206.265
+else:
+	pixel_resolution = 0
+if args.nomove:
+    nomove = True
+else:
+	nomove = False
+if args.pixelX and args.pixelY:
+    pixelX = args.pixelX
+    pixelY = args.pixelY
+    if pixel_resolution != 0:
+        radius = (max(pixelX, pixelY) * pixel_resolution) / 3600 * 1.1 #Calculate FOV and double for range of plate solving
+else:
+    pixelX = 0
+    pixelY = 0
+
+startTime = time.time()
+solveTimeFinish = time.time()
+
+#Backlash correction in arcminutes (approximately 16 full steps of each axis)
+backlashCorrection = [4.2665,0.1647]
+
+#Connect to indi server
+indiclient, blobEvent = indi.indiserverConnect()
+
+if not nomove:
+    #Initialize steppers to a known backlash position
+    initSteppers(indiclient, telescope, serialport, backlashCorrection)
+
+#Execute polar alignment routine to capture and solve three times 30 degrees apart
+p1, p1Time = slewCaptureSolve(indiclient, (targetRA % 360), targetDEC, radius, exposure, telescope, ccd, blobEvent, pixel_resolution, "capture1.fits")
+p2, p2Time = slewCaptureSolve(indiclient, ((targetRA+30) % 360), targetDEC, radius, exposure, telescope, ccd, blobEvent, pixel_resolution, "capture2.fits")
+p3, p3Time = slewCaptureSolve(indiclient, ((targetRA+60) % 360), targetDEC, radius, exposure, telescope, ccd, blobEvent, pixel_resolution, "capture3.fits")
+
+#Calculate capture time based on the average timestamps of each image
+observing_time = datetime.utcfromtimestamp(mean((p1Time, p2Time, p3Time)))
+print (f"Time of captures is {observing_time} (UTC).")
+
+#Calculate polar alignment error
+result = platesolve.polarCalc(mylat, mylong, myelev, observing_time, p1, p2, p3)
+print(f"Azimuth error correction is: {result[0]:.4f} arcminutes.")
+print(f"Altitude error correction is: {result[1]:.4f} arcminutes.")
+
+if not nomove:
+    #Disconnect telescope mount from INDI to free up serial port for Alt/Az adjustments
+    print (f"Disconnecting {telescope} from INDI server")
+    indi.disconnectScope(indiclient, telescope)
+    print ("Disconnected.")
+
+    #Adjust alt/az axis
+    print ("Adjusting altitude/azimuth axes.")
+    platesolve.adjustAltAz(result, serialport, backlashCorrection)
+
+    #Re-connect telescope mount to INDI before disconnecting from the INDI server
+    indi.connectScope(indiclient, telescope)
+
+#Disconnect from indi server
+indi.indiserverDisconnect(indiclient)
+
+print (f"Polar alignment took {time.strftime('%Mm%Ss', time.gmtime(time.time() - startTime))}.")