""" GUI tool for making initial plate estimations and manually fitting astrometric plates. """ # The MIT License # Copyright (c) 2017 Denis Vida # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal # in the Software without restriction, including without limitation the rights # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell # copies of the Software, and to permit persons to whom the Software is # furnished to do so, subject to the following conditions: # The above copyright notice and this permission notice shall be included in # all copies or substantial portions of the Software. # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN # THE SOFTWARE. from __future__ import print_function, division, absolute_import import os import sys import math import copy import time import datetime import argparse import traceback # tkinter import that works on both Python 2 and 3 try: import tkinter from tkinter import messagebox except: import Tkinter as tkinter import tkMessageBox as messagebox import numpy as np import matplotlib import matplotlib.pyplot as plt from matplotlib.font_manager import FontProperties from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes import PIL.Image, PIL.ImageDraw import scipy.optimize import scipy.ndimage from RMS.Astrometry.ApplyAstrometry import altAzToRADec, xyToRaDecPP, raDec2AltAz, raDecToXY, raDecToXYPP, \ rotationWrtHorizon, rotationWrtHorizonToPosAngle, computeFOVSize, photomLine, photometryFit, \ rotationWrtStandard, rotationWrtStandardToPosAngle, correctVignetting from RMS.Astrometry.AstrometryNetNova import novaAstrometryNetSolve from RMS.Astrometry.Conversions import date2JD, jd2Date, JD2HourAngle from RMS.Astrometry.FFTalign import alignPlatepar import RMS.ConfigReader as cr import RMS.Formats.CALSTARS as CALSTARS from RMS.Formats.Platepar import Platepar from RMS.Formats.FrameInterface import detectInputType from RMS.Formats import StarCatalog from RMS.Pickling import loadPickle, savePickle from RMS.Routines import Image from RMS.Math import angularSeparation from RMS.Misc import decimalDegreesToSexHours, openFileDialog # Import Cython functions import pyximport pyximport.install(setup_args={'include_dirs':[np.get_include()]}) from RMS.Astrometry.CyFunctions import subsetCatalog # TkAgg has issues when opening an external file prompt, so other backends are forced if available if matplotlib.get_backend() == 'TkAgg': backends = ['Qt5Agg', 'Qt4Agg'] for bk in backends: # Try setting backend try: plt.switch_backend(bk) except: pass print('Using backend: ', matplotlib.get_backend()) class FOVinputDialog(object): """ Dialog for inputting FOV centre in Alt/Az. """ def __init__(self, parent): self.parent = parent # Set initial angle values self.azim = self.alt = self.rot = 0 self.top = tkinter.Toplevel(parent) # Bind the Enter key to run the verify function self.top.bind('<Return>', self.verify) tkinter.Label(self.top, text="FOV centre (degrees) \nAzim +E of due N\nRotation from vertical").grid(row=0, columnspan=2) azim_label = tkinter.Label(self.top, text='Azim = ') azim_label.grid(row=1, column=0) self.azimuth = tkinter.Entry(self.top) self.azimuth.grid(row=1, column=1) self.azimuth.focus_set() elev_label = tkinter.Label(self.top, text='Alt =') elev_label.grid(row=2, column=0) self.altitude = tkinter.Entry(self.top) self.altitude.grid(row=2, column=1) rot_label = tkinter.Label(self.top, text='Rotation =') rot_label.grid(row=3, column=0) self.rotation = tkinter.Entry(self.top) self.rotation.grid(row=3, column=1) self.rotation.insert(0, '0') b = tkinter.Button(self.top, text="OK", command=self.verify) b.grid(row=4, columnspan=2) def verify(self, event=None): """ Check that the azimuth and altitude are withing the bounds. """ try: # Read values self.azim = float(self.azimuth.get())%360 self.alt = float(self.altitude.get()) self.rot = float(self.rotation.get())%360 # Check that the values are within the bounds if (self.alt < 0) or (self.alt > 90): messagebox.showerror(title='Range error', message='The altitude is not within the limits!') else: self.top.destroy() except: messagebox.showerror(title='Range error', message='Please enter floating point numbers, not text!') def getAltAz(self): """ Returns inputed FOV centre. """ return self.azim, self.alt, self.rot class PlateTool(object): def __init__(self, dir_path, config, beginning_time=None, fps=None, gamma=None): """ SkyFit interactive window. Arguments: dir_path: [str] Absolute path to the directory containing image files. config: [Config struct] Keyword arguments: beginning_time: [datetime] Datetime of the video beginning. Optional, only can be given for video input formats. fps: [float] Frames per second, used only when images in a folder are used. gamma: [float] Camera gamma. None by default, then it will be used from the platepar file or config. """ self.config = config self.dir_path = dir_path # If camera gamma was given, change the value in config if gamma is not None: config.gamma = gamma # Flag which regulates wheter the maxpixel or the avepixel image is shown (avepixel by default) self.img_type_flag = 'avepixel' # Star picking mode self.star_pick_mode = False self.star_selection_centroid = True self.circle_aperature = None self.circle_aperature_outer = None self.star_aperature_radius = 5 self.x_centroid = self.y_centroid = None self.closest_cat_star_indx = None self.photom_deviatons_scat = [] self.catalog_stars_visible = True self.draw_calstars = True self.draw_distorsion = False self.show_key_help = 1 # List of paired image and catalog stars self.paired_stars = [] self.residuals = None # Positions of the mouse cursor self.mouse_x = 0 self.mouse_y = 0 # Position of mouse cursor when it was last pressed self.mouse_x_press = 0 self.mouse_y_press = 0 # Kwy increment self.key_increment = 1.0 # Image gamma and levels self.bit_depth = self.config.bit_depth self.img_gamma = 1.0 self.img_level_min = self.img_level_min_auto = 0 self.img_level_max = self.img_level_max_auto = 2**self.bit_depth - 1 self.img_data_raw = None self.img_data_processed = None self.adjust_levels_mode = False self.auto_levels = False # Platepar format (json or txt) self.platepar_fmt = None # Flat field self.flat_struct = None # Dark frame self.dark = None # Image coordinates of catalog stars self.catalog_x = self.catalog_y = None self.catalog_x_filtered = self.catalog_y_filtered = None self.mag_band_string = '' # Flag indicating that the first platepar fit has to be done self.first_platepar_fit = True # Toggle zoom window self.show_zoom_window = False # Position of the zoom window (NE, NW, SE, SW) self.zoom_window_pos = '' ###################################################################################################### # Detect input file type and load appropriate input plugin self.img_handle = detectInputType(self.dir_path, self.config, beginning_time=beginning_time, fps=fps) # Extract the directory path if a file was given if os.path.isfile(self.dir_path): self.dir_path, _ = os.path.split(self.dir_path) # Load catalog stars self.catalog_stars = self.loadCatalogStars(self.config.catalog_mag_limit) self.cat_lim_mag = self.config.catalog_mag_limit # Check if the catalog exists if not self.catalog_stars.any(): messagebox.showerror(title='Star catalog error', message='Star catalog from path ' \ + os.path.join(self.config.star_catalog_path, self.config.star_catalog_file) \ + 'could not be loaded!') sys.exit() else: print('Star catalog loaded!') # Find the CALSTARS file in the given folder calstars_file = None for cal_file in os.listdir(self.dir_path): if ('CALSTARS' in cal_file) and ('.txt' in cal_file): calstars_file = cal_file break if calstars_file is None: # Check if the calstars file is required if self.img_handle.require_calstars: messagebox.showinfo(title='CALSTARS error', \ message='CALSTARS file could not be found in the given directory!') self.calstars = {} else: # Load the calstars file calstars_list = CALSTARS.readCALSTARS(self.dir_path, calstars_file) # Convert the list to a dictionary self.calstars = {ff_file: star_data for ff_file, star_data in calstars_list} print('CALSTARS file: ' + calstars_file + ' loaded!') # Load the platepar file self.platepar_file, self.platepar = self.loadPlatepar() if self.platepar_file: print('Platepar loaded:', self.platepar_file) # Print the field of view size print("FOV: {:.2f} x {:.2f} deg".format(*computeFOVSize(self.platepar))) # If the platepar file was not loaded, set initial values from config else: self.makeNewPlatepar(update_image=False) # Create the name of the platepar file self.platepar_file = os.path.join(self.dir_path, self.config.platepar_name) # Set the given gamma value to platepar if gamma is not None: self.platepar.gamma = gamma ### INIT IMAGE ### self.fig, self.ax = plt.subplots(facecolor='black') # Init the first image self.updateImage(first_update=True) # Register keys with matplotlib self.registerEventHandling() def registerEventHandling(self): """ Register mouse button and key pressess with matplotlib. """ self.fig.canvas.set_window_title('SkyFit') # Set the bacground color to black #matplotlib.rcParams['axes.facecolor'] = 'k' # Disable standard matplotlib keyboard shortcuts plt.rcParams['keymap.save'] = '' plt.rcParams['keymap.fullscreen'] = '' plt.rcParams['keymap.all_axes'] = '' plt.rcParams['keymap.quit'] = '' plt.rcParams['keymap.pan'] = '' plt.rcParams['keymap.forward'] = '' plt.rcParams['keymap.back'] = '' self.ax.figure.canvas.mpl_connect('button_press_event', self.onMousePress) self.ax.figure.canvas.mpl_connect('button_release_event', self.onMouseRelease) self.ax.figure.canvas.mpl_connect('motion_notify_event', self.onMouseMotion) # Register which mouse/keyboard events will evoke which function self.ax.figure.canvas.mpl_connect('scroll_event', self.onScroll) self.scroll_counter = 0 self.ax.figure.canvas.mpl_connect('key_press_event', self.onKeyPress) def saveState(self): """ Save the current state of the program to a file, so it can be reloaded. """ # state_date_str = datetime.datetime.now().strftime("%Y%m%d_%H%M%S.%f")[:-3] # state_file = 'skyFit_{:s}.state'.format(state_date_str) # # Save the state to a pickle file # savePickle(self, self.dir_path, state_file) # Write the latest pickle fine savePickle(self, self.dir_path, 'skyFit_latest.state') print("Saved state to file") def onMousePress(self, event): """ Called when the mouse click is pressed. """ # Record the mouse cursor positions self.mouse_x_press = event.xdata self.mouse_y_press = event.ydata def onMouseRelease(self, event): """ Called when the mouse click is released. """ # Do nothing if some button on the toolbar is active (e.g. zoom, move) if plt.get_current_fig_manager().toolbar.mode: return None # Call the same function for mouse movements to update the variables in the background self.onMouseMotion(event) # If the histogram is on, adjust the levels if self.adjust_levels_mode: # Left mouse button sets the minimum level if event.button == 1: # Compute the new image level from clicked position img_level_min_new = event.xdata/self.img_data_raw.shape[1]*(2**self.bit_depth) # Make sure the minimum level is smaller than the maximum level if (img_level_min_new < self.img_level_max) and (img_level_min_new > 0): self.img_level_min = img_level_min_new # Right mouse button sets the maximum level if event.button == 3: # Compute the new image level from clicked position img_level_max_new = event.xdata/self.img_data_raw.shape[1]*(2**self.bit_depth) # Make sure the minimum level is smaller than the maximum level if (img_level_max_new > self.img_level_min) and (img_level_max_new < (2**self.bit_depth - 1)): self.img_level_max = img_level_max_new self.updateImage() # If the star picking mode is on elif self.star_pick_mode: # Left mouse button, select stars if event.button == 1: # If the centroid of the star has to be picked if self.star_selection_centroid: # If CTRL is pressed, place the pick manually - NOTE: the intensity might be off then!!! if event.key == 'control': self.x_centroid = self.mouse_x_press self.y_centroid = self.mouse_y_press # Compute the star intensity _, _, self.star_intensity = self.centroidStar(prev_x_cent=self.x_centroid, \ prev_y_cent=self.y_centroid) # Centroid the star around the pick else: # Perform centroiding with 2 iterations x_cent_tmp, y_cent_tmp, _ = self.centroidStar() # Check that the centroiding was successful if x_cent_tmp is not None: # Centroid the star around the pressed coordinates self.x_centroid, self.y_centroid, \ self.star_intensity = self.centroidStar(prev_x_cent=x_cent_tmp, \ prev_y_cent=y_cent_tmp) else: return None # Draw the centroid on the image self.ax.scatter(self.x_centroid, self.y_centroid, marker='+', c='y', s=100, lw=3, \ alpha=0.5) # Select the closest catalog star to the centroid as the first guess self.closest_cat_star_indx = self.findClosestCatalogStarIndex(self.x_centroid, \ self.y_centroid) # Plot the closest star as a purple cross self.selected_cat_star_scatter = self.ax.scatter(self.catalog_x[self.closest_cat_star_indx], self.catalog_y[self.closest_cat_star_indx], marker='+', c='purple', s=100, lw=3) # Update canvas self.fig.canvas.draw() # Switch to the mode where the catalog star is selected self.star_selection_centroid = False self.drawCursorCircle() # If the catalog star has to be picked else: # Select the closest catalog star self.closest_cat_star_indx = self.findClosestCatalogStarIndex(self.mouse_x_press, \ self.mouse_y_press) if self.selected_cat_star_scatter is not None: self.selected_cat_star_scatter.remove() # Plot the closest star as a purple cross self.selected_cat_star_scatter = self.ax.scatter(self.catalog_x[self.closest_cat_star_indx], \ self.catalog_y[self.closest_cat_star_indx], marker='+', c='purple', s=50, lw=2) # Update canvas self.fig.canvas.draw() # Right mouse button, deselect stars elif event.button == 3: if self.star_selection_centroid: # Find the closest picked star picked_indx = self.findClosestPickedStarIndex(self.mouse_x_press, self.mouse_y_press) if self.paired_stars: # Remove the picked star from the list self.paired_stars.pop(picked_indx) self.updateImage() def onMouseMotion(self, event): """ Called with the mouse is moved. """ # Read mouse position self.mouse_x = event.xdata self.mouse_y = event.ydata # Change the position of the star aperature circle if self.star_pick_mode: self.drawCursorCircle() def drawCursorCircle(self): """ Adds a circle around the mouse cursor. """ # Delete the old aperature circle if self.circle_aperature is not None: self.circle_aperature.remove() self.circle_aperature = None if self.circle_aperature_outer is not None: self.circle_aperature_outer.remove() self.circle_aperature_outer = None # If centroid selection is on, show the annulus around the cursor if self.star_selection_centroid: # Plot a circle of the given radius around the cursor self.circle_aperature = matplotlib.patches.Circle((self.mouse_x, self.mouse_y), self.star_aperature_radius, edgecolor='yellow', fc='none') # Plot a circle of the given radius around the cursor, which is used to determine the region where the # background will be taken from self.circle_aperature_outer = matplotlib.patches.Circle((self.mouse_x, self.mouse_y), self.star_aperature_radius*2, edgecolor='yellow', fc='none', linestyle='dotted') self.ax.add_patch(self.circle_aperature) self.ax.add_patch(self.circle_aperature_outer) # If the catalog star selection mode is on, show a purple circle else: # Plot a purple circle self.circle_aperature = matplotlib.patches.Circle((self.mouse_x, self.mouse_y), 10, edgecolor='purple', fc='none') self.ax.add_patch(self.circle_aperature) # Draw the zoom window if self.show_zoom_window: self.drawZoomWindow() self.fig.canvas.draw() def drawZoomWindow(self): """ Draw the zoom window in the corner. """ if self.star_pick_mode: # If the mouse is outside the image, don't update anything if (self.mouse_x is None) or (self.mouse_y is None): return None img_w_half = self.current_ff.ncols//2 img_h_half = self.current_ff.nrows//2 # # Update the location of the zoom window # anchor_location = '' # if self.mouse_y > img_h_half: # anchor_location += 'N' # else: # anchor_location += 'S' # if self.mouse_x > img_w_half: # anchor_location += 'W' # else: # anchor_location += 'E' # self.zoom_axis.set_anchor(anchor_location) ### Extract a portion of image around the aperture ### window_radius = int(2*self.star_aperature_radius + np.sqrt(self.star_aperature_radius)) x_min_orig = int(round(self.mouse_x - window_radius)) if x_min_orig < 0: x_min = 0 else: x_min = x_min_orig x_max_orig = int(round(self.mouse_x + window_radius)) if x_max_orig > self.current_ff.ncols - 1: x_max = int(self.current_ff.ncols - 1) else: x_max = x_max_orig y_min_orig = int(round(self.mouse_y - window_radius)) if y_min_orig < 0: y_min = 0 else: y_min = y_min_orig y_max_orig = int(round(self.mouse_y + window_radius)) if y_max_orig > self.current_ff.nrows - 1: y_max = int(self.current_ff.nrows - 1) else: y_max = y_max_orig # Crop the image img_crop = np.zeros((2*window_radius, 2*window_radius), dtype=self.img_data_processed.dtype) dx_min = x_min - x_min_orig dx_max = x_max - x_max_orig + 2*window_radius dy_min = y_min - y_min_orig dy_max = y_max - y_max_orig + 2*window_radius img_crop[dy_min:dy_max, dx_min:dx_max] = self.img_data_processed[y_min:y_max, x_min:x_max] ### ### # Zoom the image zoom_factor = 8 # Make sure that the zoomed image is every larger than 1/2 of the whole image if 2*zoom_factor*window_radius > np.min([self.fig.bbox.ymax, self.fig.bbox.xmax])/2: # Compute a new zoom factor zoom_factor = np.floor((np.min([self.fig.bbox.ymax, \ self.fig.bbox.xmax])/2)/(2*window_radius)) # Don't apply zoom if the image will be smaller if zoom_factor <= 1: return None img_crop = scipy.ndimage.zoom(img_crop, zoom_factor, order=0) # Compute where the zoom will be shown on the image zoom_window_pos = '' if self.mouse_y < img_h_half: yo = 0 zoom_window_pos += 'N' else: yo = self.fig.bbox.ymax - zoom_factor*2*window_radius zoom_window_pos += 'S' if self.mouse_x > img_w_half: xo = 0 zoom_window_pos += 'W' else: xo = self.fig.bbox.xmax - zoom_factor*2*window_radius zoom_window_pos += 'E' # If the position of the zoom window has changed, reset the image if self.zoom_window_pos != zoom_window_pos: self.updateImage() self.zoom_window_pos = zoom_window_pos # Draw aperture circle to the image img = PIL.Image.fromarray(img_crop) img = img.convert('RGB') draw = PIL.ImageDraw.Draw(img) ul = zoom_factor*(window_radius - self.star_aperature_radius) br = zoom_factor*(window_radius + self.star_aperature_radius) draw.ellipse((ul, ul, br, br), outline='yellow') ul = zoom_factor*(window_radius - 2*self.star_aperature_radius) br = zoom_factor*(window_radius + 2*self.star_aperature_radius) draw.ellipse((ul, ul, br, br), outline='yellow') ### Draw centroids in the zoomed image ### def drawMarkersOnZoom(draw, x, y, color='blue', marker='x', width=1): # Check if the picked star is in the window, and plot it if it is if (x >= x_min_orig) and (x <= x_max_orig) and (y >= y_min_orig) and (y <= y_max_orig): xp = zoom_factor*(x - x_min_orig) + 1 yp = zoom_factor*(y - y_min_orig) + 1 if marker == '+': # Draw an + draw.line((xp + zoom_factor, yp, xp - zoom_factor, yp), fill=color, width=width) draw.line((xp, yp - zoom_factor, xp, yp + zoom_factor), fill=color, width=width) else: # Draw an X draw.line((xp + zoom_factor, yp + zoom_factor, xp - zoom_factor, yp - zoom_factor), \ fill=color, width=width) draw.line((xp + zoom_factor, yp - zoom_factor, xp - zoom_factor, yp + zoom_factor), \ fill=color, width=width) # Plot centroid if self.x_centroid is not None: drawMarkersOnZoom(draw, self.x_centroid, self.y_centroid, color='green', marker='x', width=2) # Plot paired stars for paired_star in self.paired_stars: img_star, catalog_star = paired_star star_x, star_y, px_intens = img_star drawMarkersOnZoom(draw, star_x, star_y, color='blue', marker='x', width=2) # Draw catalog stars if self.catalog_stars_visible: for cat_x, cat_y in zip(self.catalog_x_filtered, self.catalog_y_filtered): drawMarkersOnZoom(draw, cat_x, cat_y, color='red', marker='+') # Plot paired catalog star if self.closest_cat_star_indx is not None: drawMarkersOnZoom(draw, self.catalog_x[self.closest_cat_star_indx], \ self.catalog_y[self.closest_cat_star_indx], color='purple', marker='x', width=2) ### # Convert the PIL image object back to array img_crop = np.array(img) # Plot the zoomed image self.fig.figimage(img_crop, xo=xo, yo=yo, zorder=5, cmap='gray', \ vmin=np.min(self.img_data_processed), vmax=np.max(self.img_data_processed)) def checkParamRange(self): """ Checks that the astrometry parameters are within the allowed range. """ # Right ascension should be within 0-360 self.platepar.RA_d = self.platepar.RA_d%360 # Keep the declination in the allowed range if self.platepar.dec_d >= 90: self.platepar.dec_d = 89.999 if self.platepar.dec_d <= -90: self.platepar.dec_d = -89.999 def photometry(self, show_plot=False): """ Perform the photometry on selectes stars. """ if self.star_pick_mode: ### Make a photometry plot # Extract star intensities and star magnitudes star_coords = [] radius_list = [] px_intens_list = [] catalog_mags = [] for paired_star in self.paired_stars: img_star, catalog_star = paired_star star_x, star_y, px_intens = img_star _, _, star_mag = catalog_star # Skip intensities which were not properly calculated lsp = np.log10(px_intens) if np.isnan(lsp) or np.isinf(lsp): continue star_coords.append([star_x, star_y]) radius_list.append(np.hypot(star_x - self.platepar.X_res/2, star_y - self.platepar.Y_res/2)) px_intens_list.append(px_intens) catalog_mags.append(star_mag) # Make sure there are more than 3 stars picked if len(px_intens_list) > 3: # Fit the photometric offset (disable vignetting fit if a flat is used) photom_params, fit_stddev, fit_resids = photometryFit(px_intens_list, radius_list, \ catalog_mags, fixed_vignetting=(0.0 if self.flat_struct is not None else None)) photom_offset, vignetting_coeff = photom_params # Set photometry parameters self.platepar.mag_0 = -2.5 self.platepar.mag_lev = photom_offset self.platepar.mag_lev_stddev = fit_stddev self.platepar.vignetting_coeff = vignetting_coeff # Remove previous photometry deviation labels if len(self.photom_deviatons_scat) > 0: for entry in self.photom_deviatons_scat: resid_lbl, mag_lbl = entry try: resid_lbl.remove() mag_lbl.remove() except: pass self.photom_deviatons_scat = [] if self.catalog_stars_visible: # Plot photometry deviations on the main plot as colour coded rings star_coords = np.array(star_coords) star_coords_x, star_coords_y = star_coords.T for star_x, star_y, fit_diff, star_mag in zip (star_coords_x, star_coords_y, fit_resids, \ catalog_mags): photom_resid_txt = "{:.2f}".format(fit_diff) # Determine the size of the residual text, larger the residual, larger the text photom_resid_size = 8 + np.abs(fit_diff)/(np.max(np.abs(fit_resids))/5.0) # Plot the residual as text under the star photom_resid_lbl = self.ax.text(star_x, star_y + 10, photom_resid_txt, \ verticalalignment='top', horizontalalignment='center', \ fontsize=photom_resid_size, color='w') # Plot the star magnitude star_mag_lbl = self.ax.text(star_x, star_y - 10, "{:+6.2f}".format(star_mag), \ verticalalignment='bottom', horizontalalignment='center', \ fontsize=10, color='r') self.photom_deviatons_scat.append([photom_resid_lbl, star_mag_lbl]) self.fig.canvas.draw_idle() # Show the photometry fit plot if show_plot: ### PLOT PHOTOMETRY FIT ### # Note: An almost identical code exists in Utils.CalibrationReport # Init plot for photometry fig_p, (ax_p, ax_r) = plt.subplots(nrows=2, facecolor=None, figsize=(6.4, 7.2), \ gridspec_kw={'height_ratios':[2, 1]}) # Set photometry window title fig_p.canvas.set_window_title('Photometry') # Plot catalog magnitude vs. raw logsum of pixel intensities lsp_arr = np.log10(np.array(px_intens_list)) ax_p.scatter(-2.5*lsp_arr, catalog_mags, s=5, c='r', zorder=3, alpha=0.5, \ label="Raw") # Plot catalog magnitude vs. raw logsum of pixel intensities (only when no flat is used) if self.flat_struct is None: lsp_corr_arr = np.log10(correctVignetting(np.array(px_intens_list), \ np.array(radius_list), self.platepar.vignetting_coeff)) ax_p.scatter(-2.5*lsp_corr_arr, catalog_mags, s=5, c='b', zorder=3, alpha=0.5, \ label="Corrected for vignetting") x_min, x_max = ax_p.get_xlim() y_min, y_max = ax_p.get_ylim() x_min_w = x_min - 3 x_max_w = x_max + 3 y_min_w = y_min - 3 y_max_w = y_max + 3 # Plot fit info fit_info = "Fit: {:+.1f}*LSP + {:.2f} +/- {:.2f} ".format(self.platepar.mag_0, \ self.platepar.mag_lev, fit_stddev) \ + "\nVignetting coeff = {:.5f}".format(self.platepar.vignetting_coeff) \ + "\nGamma = {:.2f}".format(self.platepar.gamma) print(fit_info) # Plot the line fit logsum_arr = np.linspace(x_min_w, x_max_w, 10) ax_p.plot(logsum_arr, photomLine((10**(logsum_arr/(-2.5)), np.zeros_like(logsum_arr)), \ photom_offset, self.platepar.vignetting_coeff), label=fit_info, \ linestyle='--', color='k', alpha=0.5, zorder=3) ax_p.legend() ax_p.set_ylabel("Catalog magnitude ({:s})".format(self.mag_band_string)) ax_p.set_xlabel("Uncalibrated magnitude") # Set wider axis limits ax_p.set_xlim(x_min_w, x_max_w) ax_p.set_ylim(y_min_w, y_max_w) ax_p.invert_yaxis() ax_p.invert_xaxis() ax_p.grid() ### ### PLOT MAG DIFFERENCE BY RADIUS img_diagonal = np.hypot(self.platepar.X_res/2, self.platepar.Y_res/2) # Plot radius from centre vs. fit residual (including vignetting) ax_r.scatter(radius_list, fit_resids, s=5, c='b', alpha=0.5, zorder=3) # Plot a zero line ax_r.plot(np.linspace(0, img_diagonal, 10), np.zeros(10), linestyle='dashed', alpha=0.5, \ color='k') # Plot the vignetting curve (only when no flat is used) if self.flat_struct is None: # Plot radius from centre vs. fit residual (excluding vignetting fit_resids_novignetting = catalog_mags - photomLine((np.array(px_intens_list), \ np.array(radius_list)), photom_offset, 0.0) ax_r.scatter(radius_list, fit_resids_novignetting, s=5, c='r', alpha=0.5, zorder=3) px_sum_tmp = 1000 radius_arr_tmp = np.linspace(0, img_diagonal, 50) # Plot the vignetting curve vignetting_loss = 2.5*np.log10(px_sum_tmp) \ - 2.5*np.log10(correctVignetting(px_sum_tmp, radius_arr_tmp, \ self.platepar.vignetting_coeff)) ax_r.plot(radius_arr_tmp, vignetting_loss, linestyle='dotted', alpha=0.5, color='k') ax_r.grid() ax_r.set_ylabel("Fit residuals (mag)") ax_r.set_xlabel("Radius from centre (px)") ax_r.set_xlim(0, img_diagonal) fig_p.tight_layout() fig_p.show() else: print('Need more than 2 stars for photometry plot!') def updateRefRADec(self, skip_rot_update=False): """ Update the reference RA and Dec from Alt/Az. """ if not skip_rot_update: # Save the current rotation w.r.t horizon value self.platepar.rotation_from_horiz = rotationWrtHorizon(self.platepar) # Compute the datetime object of the reference Julian date time_data = [jd2Date(self.platepar.JD, dt_obj=True)] # Convert the reference alt/az to reference RA/Dec _, ra_data, dec_data = altAzToRADec(self.platepar.lat, self.platepar.lon, self.platepar.UT_corr, time_data, [self.platepar.az_centre], [self.platepar.alt_centre], dt_time=True) # Assign the computed RA/Dec to platepar self.platepar.RA_d = ra_data[0] self.platepar.dec_d = dec_data[0] if not skip_rot_update: # Update the position angle so that the rotation wrt horizon doesn't change self.platepar.pos_angle_ref = rotationWrtHorizonToPosAngle(self.platepar, \ self.platepar.rotation_from_horiz) def onKeyPress(self, event): """ Traige what happes when an individual key is pressed. """ # Switch images if event.key == 'left': self.prevImg() elif event.key == 'right': self.nextImg() elif event.key == 'm': self.toggleImageType() elif event.key == ',': # Decrement UT correction self.platepar.UT_corr -= 0.5 # Update platepar JD self.platepar.JD += 0.5/24 self.updateImage() elif event.key == '.': # Decrement UT correction self.platepar.UT_corr += 0.5 # Update platepar JD self.platepar.JD -= 0.5/24 self.updateImage() # Move RA/Dec elif event.key == 'a': self.platepar.az_centre += self.key_increment self.updateRefRADec() self.updateImage() elif event.key == 'd': self.platepar.az_centre -= self.key_increment self.updateRefRADec() self.updateImage() elif event.key == 'w': self.platepar.alt_centre -= self.key_increment self.updateRefRADec() self.updateImage() elif event.key == 's': self.platepar.alt_centre += self.key_increment self.updateRefRADec() self.updateImage() # Move rotation parameter elif event.key == 'q': self.platepar.pos_angle_ref -= self.key_increment self.updateImage() elif event.key == 'e': self.platepar.pos_angle_ref += self.key_increment self.updateImage() # Change catalog limiting magnitude elif event.key == 'r': self.cat_lim_mag += 0.1 self.catalog_stars = self.loadCatalogStars(self.cat_lim_mag) self.updateImage() elif event.key == 'f': self.cat_lim_mag -= 0.1 self.catalog_stars = self.loadCatalogStars(self.cat_lim_mag) self.updateImage() # Show/hide keyboard shortcut help elif event.key == 'f1': # Go through states of text visibility self.show_key_help += 1 if self.show_key_help >= 3: self.show_key_help = 0 self.updateImage() # Change image scale elif event.key == 'up': self.platepar.F_scale *= 1.0 + self.key_increment/100.0 self.updateImage() elif event.key == 'down': self.platepar.F_scale *= 1.0 - self.key_increment/100.0 self.updateImage() elif event.key == '1': # Increment X offset self.platepar.x_poly_rev[0] += 0.5 self.platepar.x_poly_fwd[0] += 0.5 self.updateImage() elif event.key == '2': # Decrement X offset self.platepar.x_poly_rev[0] -= 0.5 self.platepar.x_poly_fwd[0] -= 0.5 self.updateImage() elif event.key == '3': # Increment Y offset self.platepar.y_poly_rev[0] += 0.5 self.platepar.y_poly_fwd[0] += 0.5 self.updateImage() elif event.key == '4': # Decrement Y offset self.platepar.y_poly_rev[0] -= 0.5 self.platepar.y_poly_fwd[0] -= 0.5 self.updateImage() elif event.key == '5': # Decrement X 1st order distortion self.platepar.x_poly_rev[1] -= 0.01 self.platepar.x_poly_fwd[1] -= 0.01 self.updateImage() elif event.key == '6': # Increment X 1st order distortion self.platepar.x_poly_rev[1] += 0.01 self.platepar.x_poly_fwd[1] += 0.01 self.updateImage() elif event.key == '7': # Decrement Y 1st order distortion self.platepar.y_poly_rev[2] -= 0.01 self.platepar.y_poly_fwd[2] -= 0.01 self.updateImage() elif event.key == '8': # Increment Y 1st order distortion self.platepar.y_poly_rev[2] += 0.01 self.platepar.y_poly_fwd[2] += 0.01 self.updateImage() # Key increment elif event.key == '+': if self.key_increment <= 0.091: self.key_increment += 0.01 elif self.key_increment <= 0.91: self.key_increment += 0.1 else: self.key_increment += 1.0 # Don't allow the increment to be larger than 20 if self.key_increment > 20: self.key_increment = 20 self.updateImage() elif event.key == '-': if self.key_increment <= 0.11: self.key_increment -= 0.01 elif self.key_increment <= 1.11: self.key_increment -= 0.1 else: self.key_increment -= 1.0 # Don't allow the increment to be smaller than 0 if self.key_increment <= 0: self.key_increment = 0.01 self.updateImage() # Enter FOV centre elif event.key == 'v': self.platepar.RA_d, self.platepar.dec_d, self.platepar.rotation_from_horiz = self.getFOVcentre() # Recalculate reference alt/az self.platepar.az_centre, self.platepar.alt_centre = raDec2AltAz(self.platepar.JD, \ self.platepar.lon, self.platepar.lat, self.platepar.RA_d, self.platepar.dec_d) # Compute the position angle self.platepar.pos_angle_ref = rotationWrtHorizonToPosAngle(self.platepar, \ self.platepar.rotation_from_horiz) self.updateImage() # Write out the new platepar elif event.key == 'ctrl+s': # If the platepar is new, save it to the working directory if not self.platepar_file: self.platepar_file = os.path.join(self.dir_path, self.config.platepar_name) # Save the platepar file self.platepar.write(self.platepar_file, fmt=self.platepar_fmt, fov=computeFOVSize(self.platepar)) print('Platepar written to:', self.platepar_file) # Save the state self.saveState() # Save the platepar as default (SHIFT+CTRL+S) elif event.key == 'ctrl+S': platepar_default_path = os.path.join(os.getcwd(), self.config.platepar_name) # Save the platepar file self.platepar.write(platepar_default_path, fmt=self.platepar_fmt) print('Default platepar written to:', platepar_default_path) # Create a new platepar elif event.key == 'ctrl+n': self.makeNewPlatepar() # Get initial parameters from astrometry.net elif (event.key == 'ctrl+x') or (event.key == 'ctrl+X'): # Overlay text on image indicating that astrometry.net is running self.ax.text(self.img_data_raw.shape[1]/2, self.img_data_raw.shape[0]/2, \ "Solving with astrometry.net...", color='r', alpha=0.5, fontsize=16, ha='center', va='center') #self.ax.draw() self.fig.canvas.draw() self.fig.canvas.flush_events() # If shift was pressed, send only the list of x,y coords of extracted stars upload_image = True if event.key == 'ctrl+X': upload_image = False # Estimate initial parameters using astrometry.net self.getInitialParamsAstrometryNet(upload_image=upload_image) self.updateImage() # Toggle auto levels elif event.key == 'ctrl+a': self.auto_levels = not self.auto_levels self.updateImage() # Load the flat field elif event.key == 'ctrl+f': _, self.flat_struct = self.loadFlat() self.updateImage() # Load the dark frame elif event.key == 'ctrl+d': _, self.dark = self.loadDark() # Apply the dark to the flat if self.flat_struct is not None: self.flat_struct.applyDark(self.dark) self.updateImage() # Show/hide catalog stars elif event.key == 'h': self.catalog_stars_visible = not self.catalog_stars_visible self.updateImage() # Show/hide detected stars elif event.key == 'c': self.draw_calstars = not self.draw_calstars self.updateImage() # Show/hide distorsion guides elif event.key == 'ctrl+i': self.draw_distorsion = not self.draw_distorsion self.updateImage() # Increase image gamma elif event.key == 'u': # Increase image gamma by a factor of 1.1x self.updateGamma(1.1) elif event.key == 'j': # Decrease image gamma by a factor of 0.9x self.updateGamma(0.9) elif event.key == 'ctrl+h': # Toggle levels adustment mode self.adjust_levels_mode = not self.adjust_levels_mode self.updateImage() # Change modes from astrometry parameter changing to star picking elif event.key == 'ctrl+r': if self.star_pick_mode: self.disableStarPicking() self.circle_aperature = None self.circle_aperature_outer = None self.updateImage() else: self.enableStarPicking() # Toggle zoom window (SHIFT + Z) elif event.key == 'Z': if self.star_pick_mode: self.show_zoom_window = not self.show_zoom_window self.updateImage() self.drawCursorCircle() # Do a fit on the selected stars while in the star picking mode elif (event.key == 'ctrl+z') or (event.key == "ctrl+Z"): # If shift was pressed, reset distorsion parameters to zero if event.key == "ctrl+Z": self.platepar.resetDistorsionParameters() # If the first platepar is being made, do the fit twice if self.first_platepar_fit: self.fitPickedStars() self.fitPickedStars() self.first_platepar_fit = False else: # Otherwise, only fit the once self.fitPickedStars() print('Plate fitted!') elif event.key == 'enter': if self.star_pick_mode: # If the right catalog star has been selected, save the pair to the list if not self.star_selection_centroid: # Add the image/catalog pair to the list self.paired_stars.append([[self.x_centroid, self.y_centroid, self.star_intensity], self.catalog_stars[self.closest_cat_star_indx]]) # Switch back to centroiding mode self.star_selection_centroid = True self.closest_cat_star_indx = None self.updateImage() self.drawCursorCircle() elif event.key == 'escape': if self.star_pick_mode: # If the ESC is pressed when the star has been centroided, reset the centroid if not self.star_selection_centroid: self.star_selection_centroid = True self.x_centroid = None self.y_centroid = None self.star_intensity = None self.updateImage() self.drawCursorCircle() elif event.key == 'p': # Show the photometry plot self.photometry(show_plot=True) # Limit values of RA and Dec self.platepar.RA_d = self.platepar.RA_d%360 if self.platepar.dec_d > 90: self.platepar.dec_d = 90 elif self.platepar.dec_d < -90: self.platepar.dec_d = -90 def onScroll(self, event): """ Change star selector aperature on scroll. """ self.scroll_counter += event.step if self.scroll_counter > 1: self.star_aperature_radius += 1 self.scroll_counter = 0 elif self.scroll_counter < -1: self.star_aperature_radius -= 1 self.scroll_counter = 0 # Check that the star aperature is in the proper limits if self.star_aperature_radius < 2: self.star_aperature_radius = 2 if self.star_aperature_radius > 100: self.star_aperature_radius = 100 # Change the size of the star aperature circle if self.star_pick_mode: self.updateImage() self.drawCursorCircle() def toggleImageType(self): """ Toggle between the maxpixel and avepixel. """ if self.img_type_flag == 'maxpixel': self.img_type_flag = 'avepixel' else: self.img_type_flag = 'maxpixel' self.updateImage() def loadCatalogStars(self, lim_mag): """ Loads stars from the BSC star catalog. Arguments: lim_mag: [float] Limiting magnitude of catalog stars. """ # Load catalog stars catalog_stars, self.mag_band_string, self.config.star_catalog_band_ratios = StarCatalog.readStarCatalog(\ self.config.star_catalog_path, self.config.star_catalog_file, lim_mag=lim_mag, \ mag_band_ratios=self.config.star_catalog_band_ratios) return catalog_stars def drawPairedStars(self): """ Draws the stars that were picked for calibration. """ if self.star_pick_mode: # Go through all paired stars for paired_star in self.paired_stars: img_star, catalog_star = paired_star x, y, _ = img_star # Plot all paired stars self.ax.scatter(x, y, marker='x', c='b', s=100, lw=3, alpha=0.5) def drawDistorsion(self): """ Draw distorsion guides. """ # Only draw the distorsion if we have a platepar if self.platepar: # Sample 20 points on every image axis (start/end 5% from image corners) samples = 20 corner_frac = 0.05 x_samples = np.linspace(corner_frac*self.platepar.X_res, (1 - corner_frac)*self.platepar.X_res, \ samples) y_samples = np.linspace(corner_frac*self.platepar.Y_res, (1 - corner_frac)*self.platepar.Y_res, \ samples) # Create a platepar with no distorsion platepar_nodist = copy.deepcopy(self.platepar) platepar_nodist.resetDistorsionParameters() # Make X, Y pairs xx, yy = np.meshgrid(x_samples, y_samples) x_arr, y_arr = np.stack([np.ravel(xx), np.ravel(yy)], axis=-1).T # Compute RA/Dec using the normal platepar for all pairs level_data = np.ones_like(x_arr) time_data = [self.img_handle.currentTime()]*len(x_arr) _, ra_data, dec_data, _ = xyToRaDecPP(time_data, x_arr, y_arr, level_data, self.platepar) # Compute X, Y back without the distorsion jd = date2JD(*self.img_handle.currentTime()) x_nodist, y_nodist = raDecToXYPP(ra_data, dec_data, jd, platepar_nodist) # Plot the differences in X, Y data = [] color = 'r' for x0, y0, xnd, ynd in zip(x_arr, y_arr, x_nodist, y_nodist): data.append([x0, xnd]) data.append([y0, ynd]) data.append(color) plt.plot(*data, alpha=0.5) def drawLevelsAdjustmentHistogram(self, img): """ Draw a levels histogram over the image, so the levels can be adjusted. """ # If auto levels are used, show those levels if self.auto_levels: level_min = self.img_level_min_auto level_max = self.img_level_max_auto else: level_min = self.img_level_min level_max = self.img_level_max nbins = int((2**self.config.bit_depth)/2) # Compute the intensity histogram hist, bin_edges = np.histogram(img.flatten(), density=True, range=(0, 2**self.bit_depth), \ bins=nbins) # Scale the maximum histogram peak to half the image height hist *= img.shape[0]/np.max(hist)/2 # Scale the edges to image width image_to_level_scale = img.shape[1]/np.max(bin_edges) bin_edges *= image_to_level_scale ax1 = self.ax ax2 = ax1.twinx().twiny() # Plot the histogram ax2.bar(bin_edges[:-1], hist, color='white', alpha=0.5, width=img.shape[1]/nbins, edgecolor='0.5') # Plot levels limits y_range = np.linspace(0, img.shape[0], 3) x_arr = np.zeros_like(y_range) ax2.plot(x_arr + level_min*image_to_level_scale, y_range, color='w') ax2.plot(x_arr + level_max*image_to_level_scale, y_range, color='w') ax2.invert_yaxis() ax2.set_ylim([0, img.shape[0]]) ax2.set_xlim([0, img.shape[1]]) # Set background color to black ax1.set_facecolor('black') ax2.set_facecolor('black') # Plot the image level range for i in range(8): rel_i = i/8 ax2.text(rel_i*img.shape[1], 0, str(int(rel_i*2**self.bit_depth)), color='red') self.fig.sca(ax1) def updateImage(self, clear_plot=True, first_update=False): """ Update the matplotlib plot to show the current image. Keyword arguments: clear_plot: [bool] If True, the plot will be cleared before plotting again (default). first_update: [bool] Special lines will be executed if this is True, e.g. the max level will be computed. False by default. """ # Reset circle patches self.circle_aperature = None self.circle_aperature_outer = None # Limit key increment so it can be lower than 0.01 if self.key_increment < 0.01: self.key_increment = 0.01 if clear_plot: self.fig.clf() self.ax = self.fig.add_subplot(111) # Check that the calibration parameters are within the nominal range self.checkParamRange() # Load the current image self.current_ff = self.img_handle.loadChunk() if self.current_ff is None: # If the current FF couldn't be opened, go to the next messagebox.showerror(title='Read error', message='The current image is corrupted!') self.nextImg() return None # Choose appropriate image data if self.img_type_flag == 'maxpixel': img_data = self.current_ff.maxpixel else: img_data = self.current_ff.avepixel # Guess the bit depth from the array type self.bit_depth = 8*img_data.itemsize if first_update: # Set the maximum image level after reading the bit depth self.img_level_max = 2**self.bit_depth - 1 # Set the image resolution to platepar after reading the first image self.config.width = img_data.shape[1] self.config.height = img_data.shape[0] self.platepar.X_res = img_data.shape[1] self.platepar.Y_res = img_data.shape[0] # # Scale the size of the annulus (normalize so the percieved size is the same as for 1280x720) # self.star_aperature_radius *= np.sqrt(img_data.shape[1]**2 + img_data.shape[0]**2)/1500 # Apply dark if self.dark is not None: img_data = Image.applyDark(img_data, self.dark) # Apply flat if self.flat_struct is not None: img_data = Image.applyFlat(img_data, self.flat_struct) # Store image before levels modifications self.img_data_raw = np.copy(img_data) # Do auto levels if self.auto_levels: # Compute the edge percentiles self.img_level_min_auto = np.percentile(img_data, 0.1) self.img_level_max_auto = np.percentile(img_data, 99.95) # Adjust levels (auto) img_data = Image.adjustLevels(img_data, self.img_level_min_auto, self.img_gamma, \ self.img_level_max_auto, scaleto8bits=True) else: # Adjust levels (manual) img_data = Image.adjustLevels(img_data, self.img_level_min, self.img_gamma, self.img_level_max, scaleto8bits=True) # Draw levels adjustment histogram if self.adjust_levels_mode: self.drawLevelsAdjustmentHistogram(self.img_data_raw) # Store image after levels modifications self.img_data_processed = np.copy(img_data) # Show the loaded image (defining the exent speeds up image drawimg) self.ax.imshow(img_data, cmap='gray', interpolation='nearest', \ extent=(0, self.img_data_processed.shape[1], self.img_data_processed.shape[0], 0), vmin=0, vmax=255) # Draw stars that were paired in picking mode self.drawPairedStars() # Draw stars detected on this image if self.draw_calstars: self.drawCalstars() # Update centre of FOV in horizontal coordinates self.platepar.az_centre, self.platepar.alt_centre = raDec2AltAz(self.platepar.JD, self.platepar.lon, self.platepar.lat, self.platepar.RA_d, self.platepar.dec_d) ### Draw catalog stars on the image using the current platepar ### ###################################################################################################### self.catalog_x, self.catalog_y, catalog_mag = self.getCatalogStarsImagePositions(self.catalog_stars, \ self.platepar.lon, self.platepar.lat, self.platepar.RA_d, self.platepar.dec_d, \ self.platepar.pos_angle_ref, self.platepar.F_scale, self.platepar.x_poly_rev, \ self.platepar.y_poly_rev) if self.catalog_stars_visible: cat_stars = np.c_[self.catalog_x, self.catalog_y, catalog_mag] # Take only those stars inside the FOV filtered_indices, _ = self.filterCatalogStarsInsideFOV(self.catalog_stars) cat_stars = cat_stars[filtered_indices] cat_stars = cat_stars[cat_stars[:, 0] > 0] cat_stars = cat_stars[cat_stars[:, 0] < self.current_ff.ncols] cat_stars = cat_stars[cat_stars[:, 1] > 0] cat_stars = cat_stars[cat_stars[:, 1] < self.current_ff.nrows] self.catalog_x_filtered, self.catalog_y_filtered, catalog_mag_filtered = cat_stars.T if len(catalog_mag_filtered): cat_mag_faintest = np.max(catalog_mag_filtered) # Plot catalog stars self.ax.scatter(self.catalog_x_filtered, self.catalog_y_filtered, c='r', marker='+', lw=1.0, \ alpha=0.5, s=((4.0 + (cat_mag_faintest - catalog_mag_filtered))/2.0)**(2*2.512)) else: print('No catalog stars visible!') ###################################################################################################### # Draw distorsion guides if self.draw_distorsion: self.drawDistorsion() # Draw photometry if len(self.paired_stars) > 2: self.photometry() # Draw fit residuals if self.residuals is not None: self.drawFitResiduals() # Set plot limits self.ax.set_xlim([0, self.current_ff.ncols]) self.ax.set_ylim([self.current_ff.nrows, 0]) # Compute RA/Dec of the FOV centre ra_centre, dec_centre = self.computeCentreRADec() # Setup a monospace font font = FontProperties() font.set_family('monospace') font.set_size(8) if self.show_key_help == 0: text_str = 'Show fit parameters - F1' self.ax.text(10, self.current_ff.nrows, text_str, color='w', verticalalignment='bottom', \ horizontalalignment='left', fontproperties=font) # Show plate info if self.show_key_help > 0: # Show text on image with platepar parameters text_str = self.img_handle.name() + '\n' + self.img_type_flag + '\n\n' text_str += 'UT corr = {:.1f}h\n'.format(self.platepar.UT_corr) text_str += 'Ref Az = {:.3f}$\\degree$\n'.format(self.platepar.az_centre) text_str += 'Ref Alt = {:.3f}$\\degree$\n'.format(self.platepar.alt_centre) text_str += 'Rot horiz = {:.3f}$\\degree$\n'.format(rotationWrtHorizon(self.platepar)) text_str += 'Rot eq = {:.3f}$\\degree$\n'.format(rotationWrtStandard(self.platepar)) #text_str += 'Ref RA = {:.3f}\n'.format(self.platepar.RA_d) #text_str += 'Ref Dec = {:.3f}\n'.format(self.platepar.dec_d) text_str += "Scale = {:.3f}'/px\n".format(60/self.platepar.F_scale) text_str += 'Lim mag = {:.1f}\n'.format(self.cat_lim_mag) text_str += 'Increment = {:.2f}\n'.format(self.key_increment) text_str += 'Img Gamma = {:.2f}\n'.format(self.img_gamma) text_str += 'Camera Gamma = {:.2f}\n'.format(self.config.gamma) text_str += '\n' sign, hh, mm, ss = decimalDegreesToSexHours(ra_centre) if sign < 0: sign_str = '-' else: sign_str = ' ' text_str += 'RA centre = {:s}{:02d}h {:02d}m {:05.2f}s\n'.format(sign_str, hh, mm, ss) text_str += 'Dec centre = {:.3f}$\\degree$\n'.format(dec_centre) self.ax.text(10, 10, text_str, color='w', verticalalignment='top', horizontalalignment='left', \ fontproperties=font) if self.show_key_help == 1: text_str = 'Show keyboard shortcuts - F1' self.ax.text(10, self.current_ff.nrows, text_str, color='w', verticalalignment='bottom', horizontalalignment='left', fontproperties=font) # Show keyboard shortcuts if self.show_key_help > 1: # Show text on image with instructions text_str = 'Keys:\n' text_str += '-----\n' text_str += 'A/D - Azimuth\n' text_str += 'S/W - Altitude\n' text_str += 'Q/E - Position angle\n' text_str += 'Up/Down - Scale\n' text_str += '1/2 - X offset\n' text_str += '3/4 - Y offset\n' text_str += '5/6 - X 1st dist. coeff.\n' text_str += '7/8 - Y 1st dist. coeff.\n' text_str += '\n' text_str += ',/. - UT correction\n' text_str += 'R/F - Lim mag\n' text_str += '+/- - Increment\n' text_str += '\n' text_str += 'M - Toggle maxpixel/avepixel\n' text_str += 'H - Hide/show catalog stars\n' text_str += 'C - Hide/show detected stars\n' text_str += 'CTRL + I - Show/hide distorsion\n' text_str += 'U/J - Img Gamma\n' text_str += 'CTRL + H - Adjust levels\n' text_str += 'V - FOV centre\n' text_str += '\n' text_str += 'CTRL + A - Auto levels\n' text_str += 'CTRL + D - Load dark\n' text_str += 'CTRL + F - Load flat\n' text_str += 'CTRL + X - astrometry.net img upload\n' text_str += 'CTRL + SHIFT + X - astrometry.net XY\n' text_str += 'CTRL + R - Pick stars\n' text_str += 'SHIFT + Z - Show zoomed window\n' text_str += 'CTRL + N - New platepar\n' text_str += 'CTRL + S - Save platepar\n' text_str += 'SHIFT + CTRL + S - Save platepar as default\n' text_str += '\n' text_str += 'Hide on-screen text - F1\n' self.ax.text(10, self.current_ff.nrows - 5, text_str, color='w', verticalalignment='bottom', \ horizontalalignment='left', fontproperties=font) # Show fitting instructions if self.star_pick_mode: text_str = "STAR PICKING MODE" if self.show_key_help > 0: text_str += "\n'LEFT CLICK' - Centroid star\n" text_str += "'CTRL + LEFT CLICK' - Manual star position\n" text_str += "'CTRL + Z' - Fit stars\n" text_str += "'CTRL + SHIFT + Z' - Fit with initial distorsion params set to 0\n" text_str += "'P' - Photometry fit" self.ax.text(self.current_ff.ncols/2, self.current_ff.nrows, text_str, color='r', \ verticalalignment='bottom', horizontalalignment='center', fontproperties=font) self.fig.canvas.draw() def drawCalstars(self): """ Draw extracted stars on the current image. """ # Check if the given FF files is in the calstars list if self.img_handle.name() in self.calstars: # Get the stars detected on this FF file star_data = self.calstars[self.img_handle.name()] # Get star coordinates y, x, _, _ = np.array(star_data).T self.ax.scatter(x, y, edgecolors='g', marker='o', facecolors='none', alpha=0.8, \ linestyle='dotted') def updateGamma(self, gamma_adj_factor): """ Change the image gamma by a given factor. """ self.img_gamma *= gamma_adj_factor # Make sure gamma is in the proper range if self.img_gamma < 0.1: self.img_gamma = 0.1 if self.img_gamma > 10: self.img_gamma = 10 self.updateImage() def computeCentreRADec(self): """ Compute RA and Dec of the FOV centre in degrees. """ # The the time of the image img_time = self.img_handle.currentTime() # Convert the FOV centre to RA/Dec _, ra_centre, dec_centre, _ = xyToRaDecPP([img_time], [self.platepar.X_res/2], [self.platepar.Y_res/2], [1], self.platepar) ra_centre = ra_centre[0] dec_centre = dec_centre[0] return ra_centre, dec_centre def filterCatalogStarsInsideFOV(self, catalog_stars): """ Take only catalogs stars which are inside the FOV. Arguments: catalog_stars: [list] A list of (ra, dec, mag) tuples of catalog stars. """ # Get RA/Dec of the FOV centre ra_centre, dec_centre = self.computeCentreRADec() # Calculate the FOV radius in degrees fov_y, fov_x = computeFOVSize(self.platepar) fov_radius = np.sqrt(fov_x**2 + fov_y**2) # Take only those stars which are inside the FOV filtered_indices, filtered_catalog_stars = subsetCatalog(catalog_stars, ra_centre, dec_centre, \ fov_radius, self.cat_lim_mag) return filtered_indices, np.array(filtered_catalog_stars) def getCatalogStarsImagePositions(self, catalog_stars, lon, lat, ra_ref, dec_ref, pos_angle_ref, \ F_scale, x_poly_rev, y_poly_rev): """ Get image positions of catalog stars using the current platepar values. Arguments: catalog_stars: [2D list] A list of (ra, dec, mag) pairs of catalog stars. lon: [float] Longitude in degrees. lat: [float] Latitude in degrees. ra_ref: [float] Reference RA of the FOV centre (degrees). dec_ref: [float] Reference Dec of the FOV centre (degrees). pos_angle_ref: [float] Reference position angle in degrees. F_scale: [float] Image scale (px/deg). x_poly_rev: [ndarray float] Distorsion polynomial in X direction for reverse mapping. y_poly_rev: [ndarray float] Distorsion polynomail in Y direction for reverse mapping. Return: (x_array, y_array mag_catalog): [tuple of ndarrays] X, Y positons and magnitudes of stars on the image. """ ra_catalog, dec_catalog, mag_catalog = catalog_stars.T img_time = self.img_handle.currentTime() # Get the date of the middle of the FF exposure jd = date2JD(*img_time) # Convert star RA, Dec to image coordinates x_array, y_array = raDecToXY(ra_catalog, dec_catalog, jd, lat, lon, self.platepar.X_res, \ self.platepar.Y_res, ra_ref, dec_ref, self.platepar.JD, pos_angle_ref, F_scale, x_poly_rev, \ y_poly_rev, UT_corr=self.platepar.UT_corr) return x_array, y_array, mag_catalog def getPairedStarsSkyPositions(self, img_x, img_y, platepar): """ Compute RA, Dec of all paired stars on the image given the platepar. Arguments: img_x: [ndarray] Array of column values of the stars. img_y: [ndarray] Array of row values of the stars. platepar: [Platepar instance] Platepar object. Return: (ra_array, dec_array): [tuple of ndarrays] Arrays of RA and Dec of stars on the image. """ # Compute RA, Dec of image stars img_time = self.img_handle.currentTime() _, ra_array, dec_array, _ = xyToRaDecPP(len(img_x)*[img_time], img_x, img_y, len(img_x)*[1], \ platepar) return ra_array, dec_array def getInitialParamsAstrometryNet(self, upload_image=True): """ Get the estimate of the initial astrometric parameters using astromety.net. """ fail = False solution = None # Construct FOV width estimate fov_w_range = [0.75*self.config.fov_w, 1.25*self.config.fov_w] # Check if the given FF files is in the calstars list if (self.img_handle.name() in self.calstars) and (not upload_image): # Get the stars detected on this FF file star_data = self.calstars[self.img_handle.name()] # Make sure that there are at least 10 stars if len(star_data) < 10: print('Less than 10 stars on the image!') fail = True else: # Get star coordinates y_data, x_data, _, _ = np.array(star_data).T # Get astrometry.net solution, pass the FOV width estimate solution = novaAstrometryNetSolve(x_data=x_data, y_data=y_data, fov_w_range=fov_w_range) else: fail = True # Try finding the soluting by uploading the whole image if fail or upload_image: print("Uploading the whole image to astrometry.net...") # If the image is 16bit or larger, rescale and convert it to 8 bit if self.img_data_raw.itemsize*8 > 8: # Rescale the image to 8bit img_data = np.copy(self.img_data_processed) img_data -= np.min(img_data) img_data = 255*(img_data/np.max(img_data)) img_data = img_data.astype(np.uint8) else: img_data = self.img_data_raw solution = novaAstrometryNetSolve(img=img_data, fov_w_range=fov_w_range) if solution is None: messagebox.showerror(title='Astrometry.net error', \ message='Astrometry.net failed to find a solution!') return None # Extract the parameters ra, dec, orientation, scale, fov_w, fov_h = solution jd = date2JD(*self.img_handle.currentTime()) # Compute the position angle from the orientation pos_angle_ref = rotationWrtStandardToPosAngle(self.platepar, orientation) # Compute reference azimuth and altitude azim, alt = raDec2AltAz(jd, self.platepar.lon, self.platepar.lat, ra, dec) # Set parameters to platepar self.platepar.pos_angle_ref = pos_angle_ref self.platepar.F_scale = scale self.platepar.az_centre = azim self.platepar.alt_centre = alt self.updateRefRADec(skip_rot_update=True) # Save the current rotation w.r.t horizon value self.platepar.rotation_from_horiz = rotationWrtHorizon(self.platepar) # Print estimated parameters print() print('Astrometry.net solution:') print('------------------------') print(' RA = {:.2f} deg'.format(ra)) print(' Dec = {:.2f} deg'.format(dec)) print(' Azim = {:.2f} deg'.format(self.platepar.az_centre)) print(' Alt = {:.2f} deg'.format(self.platepar.alt_centre)) print(' Rot horiz = {:.2f} deg'.format(self.platepar.rotation_from_horiz)) print(' Orient eq = {:.2f} deg'.format(orientation)) print(' Pos angle = {:.2f} deg'.format(pos_angle_ref)) print(' Scale = {:.2f} arcmin/px'.format(60/self.platepar.F_scale)) def getFOVcentre(self): """ Asks the user to input the centre of the FOV in altitude and azimuth. """ # Get FOV centre root = tkinter.Tk() root.withdraw() d = FOVinputDialog(root) root.wait_window(d.top) self.azim_centre, self.alt_centre, rot_horizontal = d.getAltAz() root.destroy() # Get the middle time of the first FF img_time = self.img_handle.currentTime() # Set the reference platepar time to the time of the FF self.platepar.JD = date2JD(*img_time, UT_corr=float(self.platepar.UT_corr)) # Set the reference hour angle self.platepar.Ho = JD2HourAngle(self.platepar.JD)%360 time_data = [img_time] # Convert FOV centre to RA, Dec _, ra_data, dec_data = altAzToRADec(self.platepar.lat, self.platepar.lon, self.platepar.UT_corr, time_data, [self.azim_centre], [self.alt_centre]) return ra_data[0], dec_data[0], rot_horizontal def loadPlatepar(self): """ Open a file dialog and ask user to open the platepar file. """ platepar = Platepar() # Check if platepar exists in the folder, and set it as the default file name if it does if self.config.platepar_name in os.listdir(self.dir_path): initialfile = self.config.platepar_name else: initialfile = '' # Load the platepar file platepar_file = openFileDialog(self.dir_path, initialfile, 'Select the platepar file', matplotlib) if not platepar_file: return False, platepar # Parse the platepar file try: self.platepar_fmt = platepar.read(platepar_file, use_flat=self.config.use_flat) pp_status = True except Exception as e: print('Loading platepar failed with error:' + repr(e)) print(*traceback.format_exception(*sys.exc_info())) pp_status = False # Check if the platepar was successfuly loaded if not pp_status: messagebox.showerror(title='Platepar file error', message='The file you selected could not be loaded as a platepar file!') platepar_file, platepar = self.loadPlatepar() # Set geo location and gamma from config, if they were updated if platepar is not None: # Update the location from the config file platepar.lat = self.config.latitude platepar.lon = self.config.longitude platepar.elev = self.config.elevation # Update image resolution from config platepar.X_res = self.config.width platepar.Y_res = self.config.height # Set the camera gamma from the config file platepar.gamma = self.config.gamma # Set station ID platepar.station_code = self.config.stationID # Compute the rotation w.r.t. horizon platepar.rotation_from_horiz = rotationWrtHorizon(platepar) self.first_platepar_fit = False return platepar_file, platepar def makeNewPlatepar(self, update_image=True): """ Make a new platepar from the loaded one, but set the parameters from the config file. """ # Update the reference time img_time = self.img_handle.currentTime() self.platepar.JD = date2JD(*img_time) # Update the location from the config file self.platepar.lat = self.config.latitude self.platepar.lon = self.config.longitude self.platepar.elev = self.config.elevation # Update image resolution from config self.platepar.X_res = self.config.width self.platepar.Y_res = self.config.height # Set the camera gamma from the config file self.platepar.gamma = self.config.gamma # Estimate the scale scale_x = self.config.fov_w/self.config.width scale_y = self.config.fov_h/self.config.height self.platepar.F_scale = 1/((scale_x + scale_y)/2) # Set distorsion polynomials to zero self.platepar.x_poly_fwd *= 0 self.platepar.x_poly_rev *= 0 self.platepar.y_poly_fwd *= 0 self.platepar.y_poly_rev *= 0 # Set the first coeffs to 0.5, as that is the real centre of the FOV self.platepar.x_poly_fwd[0] = 0.5 self.platepar.x_poly_rev[0] = 0.5 self.platepar.y_poly_fwd[0] = 0.5 self.platepar.y_poly_rev[0] = 0.5 # Set station ID self.platepar.station_code = self.config.stationID # Get reference RA, Dec of the image centre self.platepar.RA_d, self.platepar.dec_d, self.platepar.rotation_from_horiz = self.getFOVcentre() # Recalculate reference alt/az self.platepar.az_centre, self.platepar.alt_centre = raDec2AltAz(self.platepar.JD, \ self.platepar.lon, self.platepar.lat, self.platepar.RA_d, self.platepar.dec_d) # Check that the calibration parameters are within the nominal range self.checkParamRange() # Compute the position angle self.platepar.pos_angle_ref = rotationWrtHorizonToPosAngle(self.platepar, \ self.platepar.rotation_from_horiz) self.platepar.auto_check_fit_refined = False # Indicate that this is the first fit of the platepar self.first_platepar_fit = True # Reset paired stars self.paired_stars = [] self.residuals = None # Indicate that a new platepar is being made self.new_platepar = True if update_image: self.updateImage() def loadFlat(self): """ Open a file dialog and ask user to load a flat field. """ # Check if flat exists in the folder, and set it as the defualt file name if it does if self.config.flat_file in os.listdir(self.dir_path): initialfile = self.config.flat_file else: initialfile = '' flat_file = openFileDialog(self.dir_path, initialfile, 'Select the flat field file', matplotlib) if not flat_file: return False, None print(flat_file) try: # Load the flat, byteswap the flat if vid file is used or UWO png flat = Image.loadFlat(*os.path.split(flat_file), dtype=self.img_data_raw.dtype, \ byteswap=self.img_handle.byteswap) except: messagebox.showerror(title='Flat field file error', \ message='Flat could not be loaded!') return False, None # Check if the size of the file matches if self.img_data_raw.shape != flat.flat_img.shape: messagebox.showerror(title='Flat field file error', \ message='The size of the flat field does not match the size of the image!') flat = None # Check if the flat field was successfuly loaded if flat is None: messagebox.showerror(title='Flat field file error', \ message='The file you selected could not be loaded as a flat field!') return flat_file, flat def loadDark(self): """ Open a file dialog and ask user to load a dark frame. """ dark_file = openFileDialog(self.dir_path, None, 'Select the dark frame file', matplotlib) if not dark_file: return False, None print(dark_file) try: # Load the dark dark = Image.loadDark(*os.path.split(dark_file), dtype=self.img_data_raw.dtype, byteswap=self.img_handle.byteswap) except: messagebox.showerror(title='Dark frame error', \ message='Dark frame could not be loaded!') return False, None dark = dark.astype(self.img_data_raw.dtype) # Check if the size of the file matches if self.img_data_raw.shape != dark.shape: messagebox.showerror(title='Dark field file error', \ message='The size of the dark frame does not match the size of the image!') dark = None # Check if the dark frame was successfuly loaded if dark is None: messagebox.showerror(title='Dark field file error', \ message='The file you selected could not be loaded as a dark field!') return dark_file, dark def nextImg(self): """ Shows the next FF file in the list. """ # Don't allow image change while in star picking mode if self.star_pick_mode: messagebox.showwarning(title='Star picking mode', message='You cannot cycle through images while in star picking mode!') return self.img_handle.nextChunk() # Reset paired stars self.paired_stars = [] self.residuals = None self.updateImage() def prevImg(self): """ Shows the previous FF file in the list. """ # Don't allow image change while in star picking mode if self.star_pick_mode: messagebox.showwarning(title='Star picking mode', message='You cannot cycle through images while in star picking mode!') return self.img_handle.prevChunk() # Reset paired stars self.paired_stars = [] self.residuals = None self.updateImage() def enableStarPicking(self): """ Enable the star picking mode where the star are manually selected for the fit. """ self.star_pick_mode = True self.updateImage() def disableStarPicking(self): """ Disable the star picking mode where the star are manually selected for the fit. """ self.star_pick_mode = False self.updateImage() def centroidStar(self, prev_x_cent=None, prev_y_cent=None): """ Find the centroid of the star clicked on the image. """ # If the centroid from the previous iteration is given, use that as the centre if (prev_x_cent is not None) and (prev_y_cent is not None): mouse_x = prev_x_cent mouse_y = prev_y_cent else: mouse_x = self.mouse_x_press mouse_y = self.mouse_y_press # Check if the mouse was pressed outside the FOV if mouse_x is None: return None, None, None ### Extract part of image around the mouse cursor ### ###################################################################################################### # Outer circle radius outer_radius = self.star_aperature_radius*2 x_min = int(round(mouse_x - outer_radius)) if x_min < 0: x_min = 0 x_max = int(round(mouse_x + outer_radius)) if x_max > self.current_ff.ncols - 1: x_max = self.current_ff.ncols - 1 y_min = int(round(mouse_y - outer_radius)) if y_min < 0: y_min = 0 y_max = int(round(mouse_y + outer_radius)) if y_max > self.current_ff.nrows - 1: y_max = self.current_ff.nrows - 1 # Crop the image img_crop = self.img_data_raw[y_min:y_max, x_min:x_max] # Perform gamma correction img_crop = Image.gammaCorrection(img_crop, self.config.gamma) ###################################################################################################### ### Estimate the background ### ###################################################################################################### bg_acc = 0 bg_counter = 0 for i in range(img_crop.shape[0]): for j in range(img_crop.shape[1]): # Calculate distance of pixel from centre of the cropped image i_rel = i - img_crop.shape[0]/2 j_rel = j - img_crop.shape[1]/2 pix_dist = math.sqrt(i_rel**2 + j_rel**2) # Take only those pixels between the inner and the outer circle if (pix_dist <= outer_radius) and (pix_dist > self.star_aperature_radius): bg_acc += img_crop[i, j] bg_counter += 1 # Calculate mean background intensity bg_intensity = bg_acc/bg_counter ###################################################################################################### ### Calculate the centroid ### ###################################################################################################### x_acc = 0 y_acc = 0 intens_acc = 0 for i in range(img_crop.shape[0]): for j in range(img_crop.shape[1]): # Calculate distance of pixel from centre of the cropped image i_rel = i - img_crop.shape[0]/2 j_rel = j - img_crop.shape[1]/2 pix_dist = math.sqrt(i_rel**2 + j_rel**2) # Take only those pixels between the inner and the outer circle if pix_dist <= self.star_aperature_radius: x_acc += j*(img_crop[i, j] - bg_intensity) y_acc += i*(img_crop[i, j] - bg_intensity) intens_acc += img_crop[i, j] - bg_intensity x_centroid = x_acc/intens_acc + x_min y_centroid = y_acc/intens_acc + y_min ###################################################################################################### return x_centroid, y_centroid, intens_acc def findClosestCatalogStarIndex(self, pos_x, pos_y): """ Finds the index of the closest catalog star on the image to the given image position. """ min_index = 0 min_dist = np.inf # Find the index of the closest catalog star to the given image coordinates for i, (x, y) in enumerate(zip(self.catalog_x, self.catalog_y)): dist = (pos_x - x)**2 + (pos_y - y)**2 if dist < min_dist: min_dist = dist min_index = i return min_index def findClosestPickedStarIndex(self, pos_x, pos_y): """ Finds the index of the closest picked star on the image to the given image position. """ min_index = 0 min_dist = np.inf picked_x = [star[0][0] for star in self.paired_stars] picked_y = [star[0][1] for star in self.paired_stars] # Find the index of the closest catalog star to the given image coordinates for i, (x, y) in enumerate(zip(picked_x, picked_y)): dist = (pos_x - x)**2 + (pos_y - y)**2 if dist < min_dist: min_dist = dist min_index = i return min_index def fitPickedStars(self): """ Fit stars that are manually picked. The function first only estimates the astrometry parameters without the distortion, then just the distortion parameters, then all together. """ # Fit the astrometry parameters, at least 5 stars are needed if len(self.paired_stars) < 4: messagebox.showwarning(title='Number of stars', message="At least 5 paired stars are needed to do the fit!") return False def _calcImageResidualsAstro(params, self, catalog_stars, img_stars): """ Calculates the differences between the stars on the image and catalog stars in image coordinates with the given astrometrical solution. """ # Extract fitting parameters ra_ref, dec_ref, pos_angle_ref, F_scale = params img_x, img_y, _ = img_stars.T # Get image coordinates of catalog stars catalog_x, catalog_y, catalog_mag = self.getCatalogStarsImagePositions(catalog_stars, \ self.platepar.lon, self.platepar.lat, ra_ref, dec_ref, pos_angle_ref, F_scale, \ self.platepar.x_poly_rev, self.platepar.y_poly_rev) # Calculate the sum of squared distances between image stars and catalog stars dist_sum = np.sum((catalog_x - img_x)**2 + (catalog_y - img_y)**2) return dist_sum def _calcSkyResidualsAstro(params, self, catalog_stars, img_stars): """ Calculates the differences between the stars on the image and catalog stars in sky coordinates with the given astrometrical solution. """ # Extract fitting parameters ra_ref, dec_ref, pos_angle_ref, F_scale = params pp_copy = copy.deepcopy(self.platepar) pp_copy.RA_d = ra_ref pp_copy.dec_d = dec_ref pp_copy.pos_angle_ref = pos_angle_ref pp_copy.F_scale = F_scale img_x, img_y, _ = img_stars.T # Get image coordinates of catalog stars ra_array, dec_array = self.getPairedStarsSkyPositions(img_x, img_y, pp_copy) ra_catalog, dec_catalog, _ = catalog_stars.T # Compute the sum of the angular separation separation_sum = np.sum(angularSeparation(np.radians(ra_array), np.radians(dec_array), \ np.radians(ra_catalog), np.radians(dec_catalog))**2) return separation_sum def _calcImageResidualsDistorsion(params, self, catalog_stars, img_stars, dimension): """ Calculates the differences between the stars on the image and catalog stars in image coordinates with the given astrometrical solution. Arguments: ... dimension: [str] 'x' for X polynomial fit, 'y' for Y polynomial fit """ if dimension == 'x': x_poly_rev = params y_poly_rev = np.zeros(12) else: x_poly_rev = np.zeros(12) y_poly_rev = params img_x, img_y, _ = img_stars.T # Get image coordinates of catalog stars catalog_x, catalog_y, catalog_mag = self.getCatalogStarsImagePositions(catalog_stars, \ self.platepar.lon, self.platepar.lat, self.platepar.RA_d, self.platepar.dec_d, \ self.platepar.pos_angle_ref, self.platepar.F_scale, x_poly_rev, y_poly_rev) # Calculate the sum of squared distances between image stars and catalog stars, per every # dimension if dimension == 'x': dist_sum = np.sum((catalog_x - img_x)**2) else: dist_sum = np.sum((catalog_y - img_y)**2) return dist_sum def _calcSkyResidualsDistorsion(params, self, catalog_stars, img_stars, dimension): """ Calculates the differences between the stars on the image and catalog stars in sky coordinates with the given astrometrical solution. Arguments: ... dimension: [str] 'x' for X polynomial fit, 'y' for Y polynomial fit """ pp_copy = copy.deepcopy(self.platepar) if dimension == 'x': pp_copy.x_poly_fwd = params else: pp_copy.y_poly_fwd = params img_x, img_y, _ = img_stars.T # Get image coordinates of catalog stars ra_array, dec_array = self.getPairedStarsSkyPositions(img_x, img_y, pp_copy) ra_catalog, dec_catalog, _ = catalog_stars.T # Compute the sum of the angular separation separation_sum = np.sum(angularSeparation(np.radians(ra_array), np.radians(dec_array), \ np.radians(ra_catalog), np.radians(dec_catalog))**2) return separation_sum # Extract paired catalog stars and image coordinates separately catalog_stars = np.array([cat_coords for img_coords, cat_coords in self.paired_stars]) img_stars = np.array([img_coords for img_coords, cat_coords in self.paired_stars]) # print('ASTRO', _calcImageResidualsAstro([self.platepar.RA_d, self.platepar.dec_d, # self.platepar.pos_angle_ref, self.platepar.F_scale], self, catalog_stars, img_stars)) # print('DIS_X', _calcImageResidualsDistorsion(self.platepar.x_poly_rev, self, catalog_stars, \ # img_stars, 'x')) # print('DIS_Y', _calcImageResidualsDistorsion(self.platepar.y_poly_rev, self, catalog_stars, \ # img_stars, 'y')) ### ASTROMETRIC PARAMETERS FIT ### # Initial parameters for the astrometric fit p0 = [self.platepar.RA_d, self.platepar.dec_d, self.platepar.pos_angle_ref, self.platepar.F_scale] # Fit the astrometric parameters using the reverse transform for reference res = scipy.optimize.minimize(_calcImageResidualsAstro, p0, args=(self, catalog_stars, img_stars), method='Nelder-Mead') # # Fit the astrometric parameters using the forward transform for reference # WARNING: USING THIS MAKES THE FIT UNSTABLE # res = scipy.optimize.minimize(_calcSkyResidualsAstro, p0, args=(self, catalog_stars, img_stars), # method='Nelder-Mead') print(res.x) # Update fitted astrometric parameters self.platepar.RA_d, self.platepar.dec_d, self.platepar.pos_angle_ref, self.platepar.F_scale = res.x # Recalculate centre self.platepar.az_centre, self.platepar.alt_centre = raDec2AltAz(self.platepar.JD, self.platepar.lon, self.platepar.lat, self.platepar.RA_d, self.platepar.dec_d) # Save the size of the image self.platepar.Y_res, self.platepar.X_res = self.current_ff.maxpixel.shape ### ### ### DISTORSION FIT ### # If there are more than 12 paired stars, fit the distortion parameters if len(self.paired_stars) > 12: ### REVERSE MAPPING FIT ### # Fit distorsion parameters in X direction, reverse mapping res = scipy.optimize.minimize(_calcImageResidualsDistorsion, self.platepar.x_poly_rev, args=(self, catalog_stars, img_stars, 'x'), method='Nelder-Mead', options={'maxiter': 10000}) # Exctact fitted X polynomial self.platepar.x_poly_rev = res.x print(res.x) # Fit distorsion parameters in Y direction, reverse mapping res = scipy.optimize.minimize(_calcImageResidualsDistorsion, self.platepar.y_poly_rev, args=(self, catalog_stars, img_stars, 'y'), method='Nelder-Mead', options={'maxiter': 10000}) # Extract fitted Y polynomial self.platepar.y_poly_rev = res.x print(res.x) ### ### # If this is the first fit of the distorsion, set the forward parametrs to be equal to the reverse if self.first_platepar_fit: self.platepar.x_poly_fwd = np.array(self.platepar.x_poly_rev) self.platepar.y_poly_fwd = np.array(self.platepar.y_poly_rev) self.first_platepar_fit = False ### FORWARD MAPPING FIT ### # Fit distorsion parameters in X direction, forward mapping res = scipy.optimize.minimize(_calcSkyResidualsDistorsion, self.platepar.x_poly_fwd, args=(self, catalog_stars, img_stars, 'x'), method='Nelder-Mead', options={'maxiter': 10000}) # Exctact fitted X polynomial self.platepar.x_poly_fwd = res.x print(res.x) # Fit distorsion parameters in Y direction, forward mapping res = scipy.optimize.minimize(_calcSkyResidualsDistorsion, self.platepar.y_poly_fwd, args=(self, catalog_stars, img_stars, 'y'), method='Nelder-Mead', options={'maxiter': 10000}) # Extract fitted Y polynomial self.platepar.y_poly_fwd = res.x print(res.x) ### ### else: print('Too few stars to fit the distorsion, only the astrometric parameters where fitted!') # Set the list of stars used for the fit to the platepar fit_star_list = [] for img_coords, cat_coords in self.paired_stars: # Compute the Julian date of the image img_time = self.img_handle.currentTime() jd = date2JD(*img_time) # Store time, image coordinate x, y, intensity, catalog ra, dec, mag fit_star_list.append([jd] + img_coords + cat_coords.tolist()) self.platepar.star_list = fit_star_list # Set the flag to indicate that the platepar was manually fitted self.auto_check_fit_refined = False ### ### ### Calculate the fit residuals for every fitted star ### # Get image coordinates of catalog stars catalog_x, catalog_y, catalog_mag = self.getCatalogStarsImagePositions(catalog_stars, \ self.platepar.lon, self.platepar.lat, self.platepar.RA_d, self.platepar.dec_d, \ self.platepar.pos_angle_ref, self.platepar.F_scale, self.platepar.x_poly_rev, \ self.platepar.y_poly_rev) residuals = [] print() print('Residuals') print('----------') print(' No, Img X, Img Y, RA (deg), Dec (deg), Mag, -2.5*LSP, Cat X, Cat Y, Err amin, Err px, Direction') # Calculate the distance and the angle between each pair of image positions and catalog predictions for star_no, (cat_x, cat_y, cat_coords, img_c) in enumerate(zip(catalog_x, catalog_y, catalog_stars, \ img_stars)): img_x, img_y, sum_intens = img_c ra, dec, mag = cat_coords delta_x = cat_x - img_x delta_y = cat_y - img_y # Compute image residual and angle of the error angle = np.arctan2(delta_y, delta_x) distance = np.sqrt(delta_x**2 + delta_y**2) # Compute the residuals in ra/dec in angular coordiniates img_time = self.img_handle.currentTime() _, ra_img, dec_img, _ = xyToRaDecPP([img_time], [img_x], [img_y], [1], self.platepar) ra_img = ra_img[0] dec_img = dec_img[0] # Compute the angular distance in degrees angular_distance = np.degrees(angularSeparation(np.radians(ra), np.radians(dec), \ np.radians(ra_img), np.radians(dec_img))) residuals.append([img_x, img_y, angle, distance, angular_distance]) # Print out the residuals print('{:3d}, {:7.2f}, {:7.2f}, {:>8.3f}, {:>+9.3f}, {:+6.2f}, {:7.2f}, {:8.2f}, {:7.2f}, {:8.2f}, {:7.2f}, {:+9.1f}'.format(star_no + 1, img_x, img_y, \ ra, dec, mag, -2.5*np.log10(sum_intens), cat_x, cat_y, 60*angular_distance, distance, np.degrees(angle))) mean_angular_error = 60*np.mean([entry[4] for entry in residuals]) # If the average angular error is larger than 60 arc minutes, report it in degrees if mean_angular_error > 60: mean_angular_error /= 60 angular_error_label = 'deg' else: angular_error_label = 'arcmin' print('Average error: {:.2f} px, {:.2f} {:s}'.format(np.mean([entry[3] for entry in residuals]), \ mean_angular_error, angular_error_label)) # Print the field of view size print("FOV: {:.2f} x {:.2f} deg".format(*computeFOVSize(self.platepar))) #################### # Save the residuals self.residuals = residuals self.updateImage() self.drawFitResiduals() def drawFitResiduals(self): """ Draw fit residuals. """ if self.residuals is not None: # Plot the residuals res_scale = 100 for entry in self.residuals: img_x, img_y, angle, distance, angular_distance = entry # Calculate coordinates of the end of the residual line res_x = img_x + res_scale*np.cos(angle)*distance res_y = img_y + res_scale*np.sin(angle)*distance # Plot the image residuals self.ax.plot([img_x, res_x], [img_y, res_y], color='orange', alpha=0.25) # Convert the angular distance from degrees to equivalent image pixels ang_dist_img = angular_distance*self.platepar.F_scale res_x = img_x + res_scale*np.cos(angle)*ang_dist_img res_y = img_y + res_scale*np.sin(angle)*ang_dist_img # Plot the sky residuals self.ax.plot([img_x, res_x], [img_y, res_y], color='yellow', alpha=0.25, linestyle='dashed') self.fig.canvas.draw_idle() if __name__ == '__main__': ### COMMAND LINE ARGUMENTS # Init the command line arguments parser arg_parser = argparse.ArgumentParser(description="Tool for fitting astrometry plates and photometric calibration.") arg_parser.add_argument('dir_path', nargs=1, metavar='DIR_PATH', type=str, \ help='Path to the folder with FF or image files, path to a video file, or to a state file. If images or videos are given, their names must be in the format: YYYYMMDD_hhmmss.uuuuuu') arg_parser.add_argument('-c', '--config', nargs=1, metavar='CONFIG_PATH', type=str, \ help="Path to a config file which will be used instead of the default one. To load the .config file in the given data directory, write '.' (dot).") arg_parser.add_argument('-t', '--timebeg', nargs=1, metavar='TIME', type=str, \ help="The beginning time of the video file in the YYYYMMDD_hhmmss.uuuuuu format.") arg_parser.add_argument('-f', '--fps', metavar='FPS', type=float, \ help="Frames per second when images are used. If not given, it will be read from the config file.") arg_parser.add_argument('-g', '--gamma', metavar='CAMERA_GAMMA', type=float, \ help="Camera gamma value. Science grade cameras have 1.0, consumer grade cameras have 0.45. Adjusting this is essential for good photometry, and doing star photometry through SkyFit can reveal the real camera gamma.") # Parse the command line arguments cml_args = arg_parser.parse_args() ######################### # If the state file was given, load the state if cml_args.dir_path[0].endswith('.state'): dir_path, state_name = os.path.split(cml_args.dir_path[0]) # Load the manual redicution tool object from a state file plate_tool = loadPickle(dir_path, state_name) # Set the dir path in case it changed plate_tool.dir_path = dir_path # Init SkyFit plate_tool.updateImage(first_update=True) plate_tool.registerEventHandling() else: # Extract the data directory path dir_path = cml_args.dir_path[0].replace('"', '') # Load the config file config = cr.loadConfigFromDirectory(cml_args.config, cml_args.dir_path) # Parse the beginning time into a datetime object if cml_args.timebeg is not None: beginning_time = datetime.datetime.strptime(cml_args.timebeg[0], "%Y%m%d_%H%M%S.%f") else: beginning_time = None # Init the plate tool instance plate_tool = PlateTool(dir_path, config, beginning_time=beginning_time, fps=cml_args.fps, gamma=cml_args.gamma) plt.tight_layout() plt.show()