Merge branch 'headers_gris_coords' into 'master'

Headers gris coords

See merge request !166

MANIFEST.in

 include kis_tools/generic/field_list/*.csv
 include kis_tools/gris/resources/*
+include kis_tools/gris/gris_coordinate_study/*.gz
 include kis_tools/util/*.pkl
 global-include *.yaml *.yml
 include kis_tools/bbi/*.csv

kis_tools/gdbs/query_gdbs.py

+import logging
+from datetime import datetime, timedelta
+
+import sqlalchemy as sa
+
+_logger = logging.getLogger(__name__)
+_gdbs_url = 'postgresql://readuser:logit@dockertest:4321/gregor'
+try:
+    gdbs_engine = sa.create_engine(_gdbs_url)
+    gdbs_engine.connect()
+except:
+    _logger.warning(
+        f'Could not connect to GREGOR database through {_gdbs_url}, are you on the KIS network and is the ssh tunnel running?')
+    raise IOError('GDBS is not online')
+
+
+def query_gdbs(table, datetime, engine=gdbs_engine):
+    query = f"""
+    SELECT value, timestamp, ttl from {table} where timestamp < '{str(datetime)}' order by timestamp desc limit 1 
+    """
+    res = engine.execute(query)
+    value, timestamp, ttl = res.fetchall()[0]
+    age = datetime - timestamp
+    is_valid = age < timedelta(seconds=ttl)
+    return value, is_valid, age
+
+
+if __name__ == '__main__':
+    engine = gdbs_engine
+    print(query_gdbs('xcen', datetime.now(), engine=engine))
+    print(query_gdbs('ycen', datetime.now(), engine=engine))

kis_tools/gris/GrisFitsFile.py

@@ -5,11 +5,11 @@ Created on 2018-04-09
 """
 
 import datetime  # date handling
+import logging
 import os  # path and file manipulations
 import sys  # error handling
 from glob import glob
 from os.path import basename, exists, join, dirname
-from warnings import warn
 
 import astropy.io.fits as fitsio  # fits reading and writing
 import astropy.units as u
@@ -18,9 +18,14 @@ import numpy as np
 from astropy.coordinates import SkyCoord, Angle
 from astropy.wcs import WCS
 
+from kis_tools.gris.coord_correction_ml import get_coords_ml
+from kis_tools.gris.util import gregor_parallactic_angle, telescope_to_hpc
+from kis_tools.util.calculations import estimate_telescope_drift
 from ..generic.fits import FitsFile
 from ..util.util import date_from_fn, gris_obs_mode
 
+logger = logging.getLogger(__name__)
+
 
 class GrisFitsFile(FitsFile):
     """Class for handling fits files for GRIS db
@@ -88,6 +93,41 @@ class GrisFitsFile(FitsFile):
         outstring += f".{int(self.runnumber):03d}"
         return outstring
 
+    @property
+    def telescope_centering(self):
+        """
+        Determine state of telescope centering at the time of observation, queries GDBS
+        for centering data
+
+        Returns:
+            outstring: x_offset, y_offset, x_uncertainty, y_uncertainty all in arcseconds
+        """
+        # Query centering from GDBS
+        try:
+            from kis_tools.gdbs.query_gdbs import query_gdbs
+            gdbs_usable = True
+            x_telcenter, x_valid, x_age = query_gdbs('xcen', self.obs_time)
+            y_telcenter, y_valid, y_age = query_gdbs('ycen', self.obs_time)
+            # Determine uncertainties for centering values
+            centering_age = max(x_age, y_age)
+            drift = estimate_telescope_drift(centering_age.seconds / 3600)
+            x_std, y_std = drift, drift
+
+            # if the centering is older than 24 hours, discard
+            if centering_age > datetime.timedelta(hours=24):
+                gdbs_usable = False
+        except IOError:
+            gdbs_usable = False
+
+        if not gdbs_usable:
+            # Determined from offset between uncentered and cross_correlated data
+            # STD-DEV of distances in X and Y again see gris coord study for details
+            uncertainty_no_center_x = 117.1
+            uncertainty_no_center_y = 182.4
+            x_telcenter, y_telcenter, x_std, y_std = 0, 0, uncertainty_no_center_x, uncertainty_no_center_y
+
+        return x_telcenter, y_telcenter, x_std, y_std
+
     @property
     def _coords_from_simple_header(self):
         """coordinates for slit, the Keywords 'SLITPOSX' and 'SLITPOSY' are assumed to be the center of the slit
@@ -104,13 +144,24 @@ class GrisFitsFile(FitsFile):
         """
 
         h = self.header
-        slit_length_pixels = h["NAXIS2"]
+        slit_length_pixels = h["NAXIS2"] if self.old_header else h['NAXIS1']
         center_pixel = slit_length_pixels // 2
 
+        # Slit center in telescope coordinates
         x_slitcenter = h["SLITPOSX"] if h["SLITPOSX"] != "" else 0
         y_slitcenter = h["SLITPOSY"] if h["SLITPOSY"] != "" else 0
-        angle = h["SLITORIE"] if h["SLITORIE"] != "" else 0
 
+        # Transform to HPC (shift by telescope centering, rotate by p0 angle)
+        delta_x, delta_y, std_delta_x, std_delta_y = self.telescope_centering
+
+        # If uncertainties are too large, use ML model to correct the coordinates
+        from kis_tools.gris.coord_correction_ml import coord_model_stdx, coord_model_stdy
+        if std_delta_x > coord_model_stdx or std_delta_y > coord_model_stdy:
+            x_hpc, y_hpc, std_delta_x, std_delta_y = get_coords_ml(self)
+        else:
+            x_hpc, y_hpc = telescope_to_hpc(x_slitcenter, y_slitcenter, self.p0_angle, delta_x, delta_y)
+
+        # Setup grid of pixels
         delta_pix = [i - center_pixel for i in range(slit_length_pixels)]
         delta_pix = np.array(delta_pix)
         arcsec_per_pixel_x = h["STEPSIZE"]  # assuming square pixels
@@ -118,19 +169,18 @@ class GrisFitsFile(FitsFile):
             h["STEPSIZE"] if h["STEPSIZE"] else 0.135
         )  # assuming square pixels
 
-        X = x_slitcenter - np.sin(angle) * delta_pix * arcsec_per_pixel_x
-        Y = y_slitcenter - np.cos(angle) * delta_pix * arcsec_per_pixel_y
+        angle = self.slit_orientation
 
-        message = (
-            "Using FITS file coordinates, they are probably wrong by at least 100 arcsec!"
-            " Also the slit angle is definitely off!"
-        )
-        warn(message)
+        # generate coordinates for full grid of slit pixels
+        X = x_hpc - np.sin(angle) * delta_pix * arcsec_per_pixel_x
+        Y = y_hpc - np.cos(angle) * delta_pix * arcsec_per_pixel_y
 
         coord_array = np.array([X, Y]) * u.arcsec
         # reorder axes for consistency with other coordinate retrieval methods
         coord_array = np.swapaxes(coord_array, 0, 1)
 
+        # track coordinate uncertainty within class
+        self.coord_uncertainty = (std_delta_x, std_delta_y)
         return coord_array
 
     @property
@@ -149,17 +199,37 @@ class GrisFitsFile(FitsFile):
             old_header = False
         return old_header
 
+    @property
+    def zenit(self):
+        """
+        Get Zenit angle (90° - elevation)
+        Returns: zenit float in degrees
+
+        """
+        zenit = 90 - self._get_value(['ELEV_ANG', 'ELEVATIO'])
+        return zenit
+
     @property
     def slit_orientation(self):
-        """Get slit_orientation, defaults to 0. Not sure what this angle means...
+        """Get slitorie:
+         clockwise angle from the western section of the solar equator
 
+        Attempts to correct for derotator and image rotation, usually correct within a couple 10s of degrees.
         Returns:
-            slit_orientation: float
+            slitorie: float (degrees)
         """
-
-        slit_orientation = self._get_value(["SLITORIE"])
-        slit_orientation = slit_orientation if slit_orientation != "" else 0
-        return slit_orientation
+        p0_angle = self.p0_angle
+        azimut = self._get_value(['AZIMUT'])
+        zenit = self.zenit
+        slitorie = self.SLITORIE
+        const_term = -130  # determined by trial and error to match results from cross correlation with HMI
+        alpha = -1 * gregor_parallactic_angle(self.obs_time) + azimut - zenit - p0_angle + const_term
+        if self.derotator_in:
+            rotangle = self.ROTANGLE
+            alpha = -2 * rotangle + (const_term / 2)
+        angle = slitorie - alpha
+
+        return angle
 
     @property
     def number_of_maps(self):
@@ -223,6 +293,15 @@ class GrisFitsFile(FitsFile):
         """
         return self._get_value(["ISTEP"])
 
+    @property
+    def p0_angle(self):
+        """
+
+        Returns:
+            p0_angle: angle between telescope and solar north
+        """
+        return self._get_value(["P0ANGLE", 'SOLAR_P0'])
+
     @property
     def step_angle(self):
         """
@@ -361,6 +440,9 @@ class GrisFitsFile(FitsFile):
         degrees = Angle(u.degree * converted[:, :2]).wrap_at(180 * u.degree)
         arcsec = degrees.to(u.arcsec)
 
+        # Track coord uncertainty within instance
+        self.coord_uncertainty = (self.CSYER1, self.CSYER2)
+
         return arcsec
 
     def plot_slit(self):

kis_tools/gris/coord_correction_ml.py

+import re
+from pathlib import Path
+
+import joblib
+import numpy as np
+
+_coord_model_path = Path(__file__).parent / 'gris_coordinate_study' / 'multi_tree_stdx38_stdy36.gz'
+_coord_model = joblib.load(_coord_model_path)
+coord_model_stdx, coord_model_stdy = (int(i) for i in
+                                      re.search(r'stdx(\d+)_stdy(\d+)', _coord_model_path.name).groups())
+
+
+def get_coords_ml(gris_fitsfile):
+    """
+    Determine the coordinates of the observation from a pre-trained statistical model.
+    The model was trained trained and evaluated against a dataset of original data and
+    results of a manual cross correlation with hmi. Albeit not being totally accurate,
+    the model has been shown to systematically reduce the distance between the given
+    and the actual coordinates. The typical accuracy (stdev) is around 40"
+    Returns: x,y tuple of the estimated observation center in Helioprojective Coordinates
+
+    """
+    features = ['SLITPOSX', 'SLITPOSY', 'AZIMUT']
+    vals = [gris_fitsfile.header[f] for f in features]
+    vals.append(gris_fitsfile.p0_angle)
+
+    # Insert bogus parallax angle, it is 0 for all cases I observed
+    # so far anyway plus is not used very strongly by the model
+    vals.insert(2, 0)
+
+    # Insert Zenit
+    vals.insert(4, gris_fitsfile.zenit)
+
+    # Insert date parameters
+    date = gris_fitsfile.obs_time
+    day_of_year = date.timetuple().tm_yday
+    vals.append(day_of_year)
+    vals.append(date.year)
+    x, y = _coord_model.predict(np.array(vals).reshape(-1, 1).T)[0]
+    return x, y, coord_model_stdx, coord_model_stdy

kis_tools/gris/gris_coordinate_study/.gitignore

-gris_coords.db

kis_tools/gris/gris_coordinate_study/__init__.py

-
-
-

kis_tools/gris/gris_coordinate_study/pull_slitorie.py

+#!/usr/bin/env python3
+
+import sys
+from pathlib import Path
+
+import pandas as pd
+
+from kis_tools.util.util import gris_run_number, date_from_fn
+
+fn = sys.argv[1]
+assert Path(fn).exists(), f'Invalid input file {fn} not found!'
+
+from scipy.io.idl import readsav
+
+data = readsav(fn)
+# rotation angle of the slit with respect to the solar equator
+angle = data['slit_angle']
+
+run = gris_run_number(fn)
+date = date_from_fn(fn)
+
+ser = pd.Series({'slit_angle': angle})
+df = pd.DataFrame()
+df[f'gris_{date.strftime("%Y%m%d")}_{run:03d}'] = ser
+df = df.T
+df.index.name = 'obs_name'
+
+from sqlalchemy import create_engine
+
+outfile = Path(__file__).parent / 'gris_coords.db'
+db_name = 'sqlite:///' + str(outfile)
+engine = create_engine(str(db_name))
+df.to_sql('slit_angle', con=engine, if_exists='append', index_label='obs_name')

kis_tools/gris/headers/wcs.py

@@ -55,6 +55,9 @@ class GrisWCSGenerator(object):
         else:
             self.coords = self.get_coords_from_header()
 
+
+
+
     def get_coords_from_header(self):
         """coordinates for slit, the Keywords 'SLITPOSX' and 'SLITPOSY' are assumed to be the center of the slit
          while the keyword 'SLITORIE' describes the rotation of the slit w.r.t. the helioprojective axes. The
@@ -66,6 +69,7 @@ class GrisWCSGenerator(object):
         """
 
         coords = self.infile.coords.value
+        self.coord_uncertainties = self.infile.coord_uncertainty
         return coords[:, 0], coords[:, 1]
 
     def get_coords_from_save(self):
@@ -193,8 +197,8 @@ class GrisWCSGenerator(object):
         # Uncertainties estimated from difference between cross correlated coordinates and coordinates
         # from the telescope. They showed a significantly larger uncertainty for the Y coordinate. See documentation of
         # https://gitlab.leibniz-kis.de/sdc/sdc_importer/issues/199
-        spatial_uncertainty_x = 5 if has_ccorr else 141
-        spatial_uncertainty_y = 5 if has_ccorr else 186
+        spatial_uncertainty_x = 5 if has_ccorr else self.coord_uncertainties[0]
+        spatial_uncertainty_y = 5 if has_ccorr else self.coord_uncertainties[1]
         c_syer = [
             spatial_uncertainty_x,
             spatial_uncertainty_y,

kis_tools/gris/study_gris_slit_orientation/clean_data.py

@@ -6,7 +6,8 @@ Created by schaffer at 11/26/19
 
 """
 import pandas as pd
-from .gris.study_gris_slit_orientation.get_data import get_raw_data
+
+from kis_tools.gris.study_gris_slit_orientation.get_data import get_raw_data
 
 
 def get_cleaned_data(*args, **kwargs):

kis_tools/gris/study_gris_slit_orientation/get_data.py

@@ -6,16 +6,17 @@ Created by schaffer at 11/26/19
 
 """
 
-from glob import iglob, glob
+from glob import glob, iglob
 from math import atan, degrees
 from os.path import abspath, dirname, basename, join
 
 from astropy.io.fits.verify import VerifyError
 from scipy.io import readsav
-from ..GrisFitsFile import GrisFitsFile
-from kis_tools.util.util import date_from_fn, gris_run_number, cache
 from tqdm import tqdm
 
+from kis_tools.util.util import date_from_fn, gris_run_number, cache
+from ..GrisFitsFile import GrisFitsFile
+
 
 def get_slit(f):
     data = readsav(f)
@@ -71,6 +72,7 @@ def angle_keywords_from_fits(save_filename):
 @cache
 def get_raw_data(force=False):
     sav_files = iglob("/dat/sdc/gris/*/context_data/*coord_mu*")
+    # sav_files = (Path("Y:") / r'dat\sdc\gris').glob("*/context_data/*coord_mu*")
     progress_bar = tqdm(list(sav_files))
     raw_values = []
 

kis_tools/gris/util.py

@@ -15,10 +15,17 @@ from itertools import chain
 from os.path import basename
 from pathlib import Path
 
+import astropy.units as u
 import numpy as np
 import pandas as pd
+import pytz
+from astroplan import Observer
+from astropy.coordinates import SkyCoord, GCRS, get_sun
 from astropy.io.fits.hdu.hdulist import fitsopen as fitsopen
+from sunpy.coordinates import Helioprojective
 
+from kis_tools.util.calculations import rotate_around_origin
+from kis_tools.util.constants import gregor_coords
 from kis_tools.util.util import gris_run_number
 
 
@@ -175,3 +182,53 @@ def extract_position_kws_from_manolo(raw_gris_file):
     )
 
     return pos
+
+
+def gregor_parallactic_angle(time, hpc_x=None, hpc_y=None):
+    """
+    Calculates parallactic angle for GREGOR based on a time and helioprojective coordinates on the sun.
+    Args:
+        time: time of observation
+        hpc_x: helioprojective x (defaults to sun center)
+        hpc_y: helioprojective y (defaults to sun center)
+
+    Returns:
+        angle: parallactic angle in degrees
+
+    """
+    # Define observer and add timezone (assume tenerife time)
+    gregor_observer = Observer(location=gregor_coords)
+    date_tz = pytz.utc.localize(time)
+    time = gregor_observer.datetime_to_astropy_time(date_tz)
+
+    # Build Sky Coord object, default to sun center at time
+    if hpc_x and hpc_y:
+        sun_coords = SkyCoord(hpc_x * u.arcsec, hpc_y * u.arcsec, frame=Helioprojective, observer='earth', obstime=time)
+        sun_coords.representation_type = 'spherical'
+        sun_coords = sun_coords.transform_to(GCRS)
+    else:
+        sun_coords = get_sun(time)
+
+    angle = gregor_observer.parallactic_angle(time, sun_coords)
+    angle = angle.to('degree').value
+    return angle
+
+
+def telescope_to_hpc(slitposx=None, slitposy=None, p0_angle=None, xcenter=None, ycenter=None):
+    """
+    Transform GRIS header coordinate data to HPC coordinates
+    Args:
+        slitposx: float SLITPOSX in arcsec
+        slitposy: float SLITPOSY in arcsec
+        p0_angle:  float p0_angle in degrees
+        xcenter: float X_center in arcsec (stored in GDBS upon centering the telescope)
+        ycenter: float Y_center in arcsec (stored in GDBS upon centering the telescope)
+
+    Returns:
+        hpc_x, hpc_y tuple of helioprojective coordinates in arcseconds
+
+    """
+    xi = slitposx - xcenter
+    yi = slitposy - ycenter
+    rotated = rotate_around_origin((xi, yi), np.radians(p0_angle))
+    return rotated

kis_tools/tests/test_coords.py

@@ -3,7 +3,6 @@
 import tempfile
 import unittest
 from itertools import combinations
-from pathlib import Path
 
 import numpy as np
 from importer_test_data import gris_structure
@@ -12,6 +11,7 @@ from scipy.io import readsav
 from kis_tools.gris.GrisDay import GrisDay
 from kis_tools.gris.GrisMapExtractor import GrisMapExtractor
 from kis_tools.gris.GrisRun import GrisRun
+from kis_tools.gris.headers.wcs import GrisWCSGenerator
 
 
 class TestCoords(unittest.TestCase):
@@ -59,3 +59,6 @@ class TestCoords(unittest.TestCase):
             if not coords_match:
                 print(msg)
             self.assertTrue(coords_match, msg=msg)
+
+    def test_gris_wcs(self):
+        wcs = GrisWCSGenerator(gris_structure.raw_files[0])

kis_tools/tests/test_gris.py

@@ -23,6 +23,7 @@ from kis_tools.gris.calib_settings.parsing import (
     get_gris_calib_scripts,
     parse_file,
 )
+from kis_tools.gris.coord_correction_ml import get_coords_ml
 from kis_tools.gris.headers.wcs import reorder_gris_data
 from kis_tools.util.util import groupby_gris_run
 
@@ -100,7 +101,6 @@ class TestGris(unittest.TestCase):
         self.assertIsInstance(gff, GrisFitsFile, "Instantiation of GrisFitsFile failed")
         # test parse function
         gff.parse()
-        gff.parse()
 
         # test printing of object
         self.assertIsInstance(gff.__str__(), str)
@@ -122,12 +122,16 @@ class TestGris(unittest.TestCase):
         gf2 = GrisFitsFile(test_data.gris.l1_data[0])
         gf2.slit_orientation
         gf2.parse()
+        gf2.telescope_centering
         _ = gf2.modestring
         _ = gf2._coords_from_simple_header
 
         # test coordinates from WCS
         self.assertEqual(len(gf2._coords_from_wcs.shape), 2)
 
+        # test coods via ml:
+        _ = get_coords_ml(gf2)
+
         # test coordinates from translated and untranslated file
         translated_file = self.filetarget
         raw_file = raw_files[0]

kis_tools/util/calculations.py

@@ -147,3 +147,21 @@ def rotate_around_origin(xy, angle_rad):
     yy = -x * np.sin(angle_rad) + y * np.cos(angle_rad)
 
     return xx, yy
+
+
+def estimate_telescope_drift(hours_since_center):
+    """
+    Gives an estimate of the drift of the telescope centering information. The parametrization of the degradation is
+    based on studies performed by Olivier Grassin and backed up by a validation by Carl Schaffer based on ~500 cross
+    correlations between GREGOR data and HMI. Further information can be found in the kis_tools repository under
+    GRIS coordinate study
+    Args:
+        hours_since_center: hours since the centering procedure was done
+
+    Returns:
+        degradation: float mean drift in arcseconds
+
+    """
+    offset, slope = 23.4, 21.7
+    degradation = offset + slope * hours_since_center
+    return degradation

requirements.txt

@@ -8,4 +8,7 @@ pytest-cov>=2.7
 scipy>=1.3
 sunpy[all]>=2
 tqdm
+joblib
+scikit-learn
 urllib3
+astroplan
\ No newline at end of file