Commit 2fbc64bc authored by vertighel's avatar vertighel
Browse files

Added custom format to avoid rounding errors

parent 2c9dd8ab
Loading
Loading
Loading
Loading
Loading
+79 −0
Original line number Diff line number Diff line
from astropy.io import fits
import contextlib

_original_astropy_format_float = None

class FormattedFloat(float):
    """
    A float subclass that carries a desired number of decimal digits
    for its string representation.
    """
    
    def __new__(cls, value, num_digits=None):
        obj = super().__new__(cls, value)
        obj._num_digits = num_digits
        return obj

    def __str__(self):
        """
        Returns the float value as a formatted string based on _num_digits.
        If _num_digits is None, falls back to default float string representation.
        """
        
        if self._num_digits is not None:
            # Use f-string formatting for fixed-point notation
            return f"{self.__float__():.{self._num_digits}f}"
        else:
            # Fallback to default float string representation
            return super().__str__()

    def __repr__(self):
        """
        Provides a developer-friendly representation of the FormattedFloat.
        """
        
        return f"FormattedFloat({super().__repr__()}, num_digits={self._num_digits})"


@contextlib.contextmanager
def custom_float():
    """
    Context manager to temporarily override Astropy's internal float formatter.
    When active, it will use the __str__ method of FormattedFloat instances
    for FITS header values.
    """
    
    global _original_astropy_format_float
    _original_astropy_format_float = fits.card._format_float

    def _new_format_float_for_fits(value):
        """
        This function replaces Astropy's default float formatter.
        If the value is a FormattedFloat, it uses its __str__ method.
        Otherwise, it delegates to the original Astropy formatter.
        """
        
        if isinstance(value, FormattedFloat):
            return str(value) # This calls FormattedFloat.__str__
        else:
            return _original_astropy_format_float(value)

    fits.card._format_float = _new_format_float_for_fits
    try:
        # Execute the code within the 'with' block        
        yield 
    finally:
        # Restore the original formatter when exiting the 'with' block
        fits.card._format_float = _original_astropy_format_float


def header_round(value, num_digits):
    """
    Returns a FormattedFloat instance that, when assigned to an Astropy FITS header
    within the 'custom_fits_float_formatter_context', will be formatted
    with the specified number of decimal digits.

    This function does NOT perform rounding itself, but rather attaches
    formatting metadata to the float for Astropy's FITS writer.
    """
    return FormattedFloat(value, num_digits)
+21 −20
Original line number Diff line number Diff line
#------------------------------------------------------------------------------|

#--------------------------------------------------------------------------------|

[Initial]
SIMPLE = True         | Standard FITS format
BITPIX = 16           | Array data type
@@ -46,6 +48,8 @@ INSTRUME = | Instrument name
OBSTYPE =             | Observation type
                      
[Detector]            
IMAGETYP =            | Frame type
FILTER =              | Photometric filter used
DETECTOR =            | Detector identifier
DETSIZE =             | [px] Physical CCD dimensions
DETROT =              | [deg] Rotation offset of the detector
@@ -56,9 +60,6 @@ SET-TEMP = | [C] CCD temperature set point
CCD-TEMP =            | [C] CCD temperature
GAIN =                | [e-/ADU] Gain
RDNOISE =             | [e- RMS] Readout noise

IMAGETYP =            | Frame type
FILTER =              | Photometric filter used
EXPTIME =             | [s] Exposure time
XBINNING =            | Binning factor in X
YBINNING =            | Binning factor in Y
+99 −0
Original line number Diff line number Diff line

#---------------------|-type--|%.f|----------------------------------------------|

[Initial]
SIMPLE = True         | bool  |   | Standard FITS format
BITPIX = 16           | int   |   | Array data type
NAXIS = 2             | int   |   | Number of data axes
NAXIS1 =              | int   |   | [px] Length of data axis 1 (X)
NAXIS2 =              | int   |   | [px] Length of data axis 2 (Y)
EXTEND = True         | bool  |   | FITS file may contain extensions
BSCALE = 1            | int   |   | Scale factor applied to data
BZERO = 32768         | int   |   | Offset applied to data after scaling
OBSERVER = UNKNOWN    | str   |   | Observer name
                                  
[Object]                          
OBJECT =              | str   |   | Name of observed object
RA =                  | str   |   | [hh:mm:ss.ss] Right Ascension, sexagesimal
DEC =                 | str   |   | [+dd:mm:ss.ss] Declination, sexagesimal
RA_DEG =              | float | 7 | [deg] Right Ascension in decimal degrees
DEC_DEG =             | float | 7 | [deg] Declination in decimal degrees
ALT =                 | float | 7 | [deg] Altitude of object above horizon
AZ =                  | float | 7 | [deg] Azimuth of object from North (E=90)
HA =                  | float | 2 | [h] Hour Angle
LST =                 | float | 2 | [h] Local Sidereal Time
PARANGLE =            | float | 2 | [deg] Object Parallactic angle
POSANGLE =            | float | 2 | [deg] Object Position angle
AIRMASS =             | float | 2 | Approximate air mass
                                  
[Time]                            
DATE =                | str   |   | [YYYY-MM-DD] File creation date
DATE-OBS =            | str   |   | [YYYY-MM-DDTHH:MM:SS] UTC observation date
MJD-OBS =             | float | 7 | [d] Modified Julian Date of observation
                                  
[Telescope]                       
TELESCOP =            | str   |   | Telescope name
FOCALLEN =            | int   |   | [mm] Telescope focal length
APTDIA =              | int   |   | [mm] Telescope aperture diameter
OBS-LONG =            | float | 7 | [deg] Observatory longitude (East > 0)
OBS-LAT =             | float | 7 | [deg] Observatory latitude (North > 0)
OBS-ELEV =            | int   |   | [m] Observatory altitude above sea level
FOCUSPOS =            | int   |   | [um] Focuser position
FOCUSTEM =            | float | 2 | [C] Focuser temperature
DEROTANG = 0          | float | 2 | [deg] Rotator angle, if any
PIERSIDE = UNKNOWN    | str   |   | Side of the camera wrt pier, if any
                                  
[Imaging]                         
INSTRUME =            | str   |   | Instrument name
OBSTYPE =             | str   |   | Observation type
                                  
[Detector]                        
IMAGETYP =            | str   |   | Frame type
FILTER =              | str   |   | Photometric filter used
DETECTOR =            | str   |   | Detector identifier
DETSIZE =             | str   |   | [px] Physical CCD dimensions
DETROT =              | float | 2 | [deg] Rotation offset of the detector
XPIXSZ =              | float | 1 | [um] Pixel X axis size
YPIXSZ =              | float | 1 | [um] Pixel Y axis size
PIXSCALE =            | float | 2 | [arcsec/px] Plate scale in binning 1
SET-TEMP =            | float | 1 | [C] CCD temperature set point
CCD-TEMP =            | float | 1 | [C] CCD temperature
GAIN =                | float | 2 | [e-/ADU] Gain
RDNOISE =             | float | 2 | [e- RMS] Readout noise
EXPTIME =             | float | 1 | [s] Exposure time
XBINNING =            | int   |   | Binning factor in X
YBINNING =            | int   |   | Binning factor in Y
                                  
[Ambient]                         
SUNALT =              | float | 1 | [deg] Sun altitude (< 0 means night)
MOONDIST =            | float | 1 | [deg] Moon-target angular distance
MOONPHAS =            | float | 2 | Lunar phase (0=new, 1=full)
TEMPERAT =            | float | 1 | [C] Ambient temperature
                                  
[WCS]                             
WCSAXES = 2           | int   |   | Number of axes
CTYPE1 = RA---TAN     | str   |   | Right ascension, gnomonic projection
CTYPE2 = DEC--TAN     | str   |   | Declination, gnomonic projection
CUNIT1 = deg          | str   |   | Units of coordinate axis 1
CUNIT2 = deg          | str   |   | Units of coordinate axis 2
CRPIX1 =              | float | 1 | [px] X reference pixel
CRPIX2 =              | float | 1 | [px] Y reference pixel
CRVAL1 =              | float | 7 | [deg] RA at reference pixel
CRVAL2 =              | float | 7 | [deg] DEC at reference pixel
CDELT1 =              | float | 7 | [deg/px] RA pixel scale
CDELT2 =              | float | 7 | [deg/px] DEC pixel scale
PC1_1 = -1.0          | float | 5 | Rotation matrix element
PC1_2 = 0.0           | float | 5 | Rotation matrix element
PC2_1 = 0.0           | float | 5 | Rotation matrix element
PC2_2 = 1.0           | float | 5 | Rotation matrix element
RADESYS = ICRS        | str   |   | Coordinate system reference frame
EQUINOX = 2000.0      | float | 1 | [yr] Equinox of coordinates
                                  
[Final]                           
FILEORIG =            | str   |   | Original file name
SWCREATE =            | str   |   | Software that created FILEORIG
CHECKSUM =            | str   |   | Checksum of header
DATASUM =             | str   |   | Data sum of FITS file
ORIGIN = NOCHE v0.1   | str   |   | Origin of FITS file
COMMENT = NOCTIS common header
HISTORY = Updated by NOCHE
+88 −55
Original line number Diff line number Diff line
@@ -13,6 +13,8 @@ from astropy.io import fits
from astropy import units as u
from astropy import log

from .custom_rounding import custom_float, header_round

class Noche:
    """
    This class provides tools to initialize and populate a FITS header for
@@ -102,6 +104,7 @@ class Noche:
            for key, value in config.items(section):
                try:
                    val, comment = value.split("|")
#                    val, typ, fmt, comment = value.split("|")
                    comment = comment.strip()
                except ValueError:
                    # HISTORY and COMMENT have no comment
@@ -349,12 +352,14 @@ class Noche:
            time = Time(obstime)
            self.set_obstime(time)

        with custom_float():
            
            self.header['RA'] = coord.ra.to_string(unit=u.hourangle, sep=':',
                                                   pad=True, precision=1)
            self.header['DEC'] = coord.dec.to_string(unit=u.deg, sep=':',
                                                     pad=True, precision=1)
        self.header['RA_DEG'] = coord.ra.deg
        self.header['DEC_DEG'] = coord.dec.deg
            self.header['RA_DEG'] = header_round(coord.ra.deg, 7)
            self.header['DEC_DEG'] = header_round(coord.dec.deg,7)

        self._update()

@@ -378,29 +383,31 @@ class Noche:
        altaz = self._coord.altaz

        # Altitudine and Azimuth
        self.header['ALT'] = round(altaz.alt.deg, 7)
        self.header['AZ'] = round(altaz.az.deg, 7)
        self.header['AIRMASS'] = round(altaz.secz.value, 2)
        with custom_float():
            
            self.header['ALT'] = header_round(altaz.alt.deg, 7)
            self.header['AZ'] = header_round(altaz.az.deg, 7)
            self.header['AIRMASS'] = header_round(altaz.secz.value, 2)

            # Local Sideral Time and Hour Angle
            lst = self._obstime.sidereal_time('mean', longitude=self._location.lon)
            ha = (lst - self._coord.ra).hour
            lst_hours = lst.hour

        self.header['LST'] = round(lst_hours, 2)
        self.header['HA'] = round(ha, 2)
            self.header['LST'] = header_round(lst_hours, 2)
            self.header['HA'] = header_round(ha, 2)

            # Position angle: with respect to Celestial North Pole
            north_celestial = SkyCoord(ra=0*u.deg, dec=90*u.deg, frame='icrs')
            posangle = self._coord.position_angle(north_celestial).to(u.deg).value
        self.header['POSANGLE'] = round(posangle, 2)
            self.header['POSANGLE'] = header_round(posangle, 2)

            # Parallactic angle:  between local meridian and celestial axis
            parangle = (posangle - altaz.az.deg + 360) % 360
            if parangle > 180:
                parangle -= 360  # Wrap to [-180, 180]
                
        self.header['PARANGLE'] = round(parangle, 2)
            self.header['PARANGLE'] = header_round(parangle, 2)


    def set_wcs(self, angle=None):
@@ -429,9 +436,10 @@ class Noche:
        _, xsize = map(int, x_str.split(':'))
        _, ysize = map(int, y_str.split(':'))
        
        with custom_float():
            crpix = [xsize/self.header["XBINNING"]/2, ysize/self.header["YBINNING"]/2]
        cdelt1 = round(self.header["PIXSCALE"]*self.header["XBINNING"]*u.arcsec.to(u.deg), 7)
        cdelt2 = round(self.header["PIXSCALE"]*self.header["YBINNING"]*u.arcsec.to(u.deg), 7)
            cdelt1 = header_round(self.header["PIXSCALE"]*self.header["XBINNING"]*u.arcsec.to(u.deg), 7)
            cdelt2 = header_round(self.header["PIXSCALE"]*self.header["YBINNING"]*u.arcsec.to(u.deg), 7)

        if not angle:
            angle = self.header["DEROTANG"] + self.header["DETROT"]
@@ -446,12 +454,13 @@ class Noche:

        flip = -1 # East to the left
        
        self.header['CRPIX1'] = crpix[0] 
        self.header['CRPIX2'] = crpix[1]
        self.header['CRVAL1'] = crval_ra
        self.header['CRVAL2'] = crval_dec
        self.header['CDELT1'] = cdelt1 * flip
        self.header['CDELT2'] = cdelt2
        with custom_float():
            self.header['CRPIX1'] = header_round(crpix[0], 7)
            self.header['CRPIX2'] = header_round(crpix[1], 7)
            self.header['CRVAL1'] = header_round(crval_ra, 7)
            self.header['CRVAL2'] = header_round(crval_dec, 7)
            self.header['CDELT1'] = header_round(cdelt1 * flip, 7)
            self.header['CDELT2'] = header_round(cdelt2, 7)
            self.header["PC1_1"] = +np.cos(angle)
            self.header["PC1_2"] = -np.sin(angle)
            self.header["PC2_1"] = +np.sin(angle)
@@ -482,18 +491,21 @@ class Noche:
            moon = get_body("moon", time)
            moon.location = loc


        with custom_float():

            # MOONDIST: angular distance between target and Moon
            moondist = moon.separation(self._coord, origin_mismatch="ignore").deg
        self.header['MOONDIST'] = round(moondist, 1)
            self.header['MOONDIST'] = header_round(moondist, 1)

            # MOONPHAS: Moon phase
            elongation = moon.separation(sun).deg
            moonphas = (1 + np.cos(np.radians(elongation))) / 2
        self.header['MOONPHAS'] = round(moonphas, 2)
            self.header['MOONPHAS'] = header_round(moonphas, 2)

            # SUNALT: Sun altitude above the horizon
            sunalt = sun.altaz.alt.deg
        self.header['SUNALT'] = round(sunalt, 1)
            self.header['SUNALT'] = header_round(sunalt, 1)


    def check_empty(self):
@@ -549,7 +561,8 @@ class Noche:
            self.set_wcs()

            
    def _parse(self, val):
            
    def _parse(self, val, typ=None, fmt=None):
        """
        Parse values from configuration file.

@@ -590,3 +603,23 @@ class Noche:
                    pass

        return val

    

    def _parse2(self, val, typ, fmt):
        type_mapping = {
            "bool": bool,
            "float": float,
            "int": int,
            "str": str,
        }
        
        if typ in type_mapping:
            the_type = type_mapping[typ.strip()]
            val = the_type(val)
            if typ.strip() == "float":
                val = header_round(val, fmt)                        
        else:
            log.error(f"Not supported type: {typ}")

        return val