From 9d5df7d74db49b308dc397fbdea5806167c414f2 Mon Sep 17 00:00:00 2001
From: Paul Chote
Date: Wed, 31 Aug 2016 17:53:23 +0100
Subject: [PATCH] Add image generation.
---
pmfinding.py | 133 ++++++++++++++++++++++++++++++++++++++++++++++++++-
1 file changed, 131 insertions(+), 2 deletions(-)
diff --git a/pmfinding.py b/pmfinding.py
index 014ea2f..f4a964d 100644
--- a/pmfinding.py
+++ b/pmfinding.py
@@ -1,6 +1,130 @@
+#
+# make-pm-findingchart is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# make-pm-findingchart is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with make-pm-findingchart. If not, see .
+
+# pylint: disable=invalid-name
+
+import argparse
+import datetime
+import io
+import math
+import urllib.request
+import numpy
+import os
+import sys
+from PIL import Image, ImageOps, ImageDraw, ImageFont
+from astropy import wcs
+from astropy.io import fits
+
from flask import Flask
from flask import render_template
from flask import request
+from flask import send_file
+
+def rescale_image_data(data, clip_low, clip_high):
+ """ Returns a normalised array where clip_low percent of the pixels are 0 and
+ clip_high percent of the pixels are 255
+ """
+ high = numpy.percentile(data, clip_high)
+ low = numpy.percentile(data, clip_low)
+ scale = 255. / (high - low)
+ data = numpy.clip(data, low, high)
+ return 255 - scale * (data - low)
+
+def parse_sexagesimal(string):
+ """Converts a sexagesimal string to decimal"""
+ parts = string.split(':')
+ if len(parts) != 3:
+ raise ValueError('Invalid input: ' + string)
+
+ a = float(parts[0])
+ b = math.copysign(float(parts[1]), a)
+ c = math.copysign(float(parts[2]), a)
+
+ return a + b / 60 + c / 3600
+
+def offset_proper_motion(ra_degrees, dec_degrees, pm_ra_degrees, pm_dec_degrees, delta_yr):
+ ra = ra_degrees + float(pm_ra_degrees) / math.cos(dec_degrees * math.pi / 180) * delta_yr
+ dec = dec_degrees + float(pm_dec_degrees) * delta_yr
+ return (ra, dec)
+
+def generate_finding_chart(out_year, in_ra, in_dec, in_year, ra_pm, dec_pm, width, height, survey):
+ circle_r = 10
+ circle_r2 = 10
+ ra_j2000_degrees = parse_sexagesimal(in_ra) * 15
+ dec_j2000_degrees = parse_sexagesimal(in_dec)
+ ra_pm_degrees = float(ra_pm) / 3600
+ dec_pm_degrees = float(dec_pm) / 3600
+
+ ra_target, dec_target = offset_proper_motion(ra_j2000_degrees, dec_j2000_degrees,
+ ra_pm_degrees, dec_pm_degrees,
+ float(out_year) - float(in_year))
+ url = 'http://archive.stsci.edu/cgi-bin/dss_search?r=' + str(ra_target) + '&dec=' + str(dec_target) \
+ + '&v=' + survey + '&f=dss1&s=on&e=J2000&h=' + str(height) + '&w=' + str(width)
+ filename, _ = urllib.request.urlretrieve(url)
+
+ try:
+ with fits.open(filename) as hdulist:
+ frame = hdulist[0]
+ arcsec_per_px = frame.header['PLTSCALE'] * frame.header['XPIXELSZ'] / 1000
+
+ # Headers can contain bogus time values (e.g. 93 minutes), so only consider year part
+ frame_date = datetime.datetime.strptime(frame.header['DATE-OBS'][0:11], '%Y-%m-%dT')
+ frame_year = (float(frame_date.strftime("%j"))-1) / 366 + float(frame_date.strftime("%Y"))
+ delta_years = frame_year - 2000
+
+ frame_coords = offset_proper_motion(ra_j2000_degrees, dec_j2000_degrees,
+ ra_pm_degrees, dec_pm_degrees,
+ float(frame_year) - float(in_year))
+
+ w = wcs.WCS(hdulist[0].header)
+ old_x, old_y = w.wcs_world2pix(numpy.array([frame_coords], numpy.float_), 0, ra_dec_order=True)[0]
+ new_x, new_y = w.wcs_world2pix(numpy.array([[ra_target, dec_target]], numpy.float_), 0, ra_dec_order=True)[0]
+
+ delta_x = new_x - old_x
+ delta_y = new_y - old_y
+ delta_l = math.sqrt(delta_x * delta_x + delta_y * delta_y)
+ dir_x = delta_x / delta_l
+ dir_y = delta_y / delta_l
+
+
+ scaled = rescale_image_data(frame.data, 1, 99.5)
+ png = Image.fromarray(scaled).convert('RGB').resize((512, 512), Image.BICUBIC)
+ scale_x = float(png.width) / scaled.shape[0]
+ scale_y = float(png.height) / scaled.shape[1]
+
+ line_start_x = old_x + circle_r2 * dir_x / scale_y
+ line_start_y = old_y + circle_r2 * dir_y / scale_y
+ line_end_x = new_x - circle_r * dir_x / scale_y
+ line_end_y = new_y - circle_r * dir_y / scale_y
+
+ arrow_a_x = line_end_x - 10 * (dir_y + dir_x) / scale_x
+ arrow_a_y = line_end_y + 10 * (-dir_y + dir_x) / scale_y
+ arrow_b_x = line_end_x - 10 * (-dir_y + dir_x) / scale_x
+ arrow_b_y = line_end_y + 10 * (-dir_y - dir_x) / scale_y
+
+ draw = ImageDraw.Draw(png)
+ draw.ellipse((scale_x * new_x-circle_r, scale_y * new_y-circle_r, scale_x * new_x + circle_r, scale_y * new_y + circle_r), fill='red')
+ draw.ellipse((scale_x * old_x-circle_r2, scale_y * old_y-circle_r2, scale_x * old_x + circle_r2, scale_y * old_y + circle_r2), outline='blue')
+
+ if delta_l * scale_x > circle_r + circle_r2:
+ draw.line((scale_x * line_start_x, scale_y * line_start_y, scale_x * line_end_x, scale_y * line_end_y), 'blue')
+ draw.line((scale_x * arrow_a_x, scale_y * arrow_a_y, scale_x * line_end_x, scale_y * line_end_y), 'blue')
+ draw.line((scale_x * arrow_b_x, scale_y * arrow_b_y, scale_x * line_end_x, scale_y * line_end_y), 'blue')
+
+ return ImageOps.flip(png)
+ finally:
+ os.remove(filename)
app = Flask(__name__)
@@ -10,5 +134,10 @@ def input_display():
@app.route('/generate')
def generate_chart():
- print(request.args.get('coordinates', ''))
- return 'hi'
+ print(request.args)
+ chart = generate_finding_chart(request.args['outepoch'], request.args['ra'], request.args['dec'], request.args['epoch'], request.args['rapm'], request.args['decpm'], request.args['size'], request.args['size'], request.args['survey'])
+ output = io.BytesIO()
+ chart.save(output, format='PNG')
+ output.seek(0)
+
+ return send_file(output, attachment_filename='chart.png', mimetype='image/png')