<commit_before>#!/usr/bin/env python3

"""
Extracts a spectrum from a region in a spectral cube.
The pixels will be averaged spatially inside the given region.
The region must be specified as shape,coords,frame,parameters
Region can be, point, box, circle or ellipse. Also, CASA region files are supported.
Coords can be pix or sky.
Frames supported include: fk4, fk5, icrs, and gal.
Parameters, for point the coordinates of the point, e.g. point,pix,,512,256.
For box the bottom left corner and top right corner coordinates, e.g. box,pix,,256,256,512,512.
For circle the center of the circle and the radius, e.g. circle,pix,,512,256,10.
For ellipse, its center, major and minor axes and the angle in degrees, e.g., 
ellipse,sky,fk5,13h05m27.462s,-49d28m06.547s,0.345293s,0.493276s,0
If the coordinates are given in sky values, then the units must included, e.g. 
point,sky,12h,2d will extract the spectrum for the pixel located at RA 12 hours and DEC 2 degrees.
For circle the radius units in sky coordinates can be either d for degrees, m for arcminutes or 
s for arcseconds. The conversion will use the scale in the RA direction of the cube, i.e. will 
use CDELT1 to convert from angular units to pixels.
"""

import os
import matplotlib as mpl
havedisplay = "DISPLAY" in os.environ
if not havedisplay:
    mpl.use('Agg')

import sys
import re
import argparse
import logging

from astropy.io import fits
from astropy.table import Table
from astropy.io import ascii
from astropy.wcs import WCS, WCSSUB_SPECTRAL
from astropy.coordinates import SkyCoord

import crrlpy.imtools as ci
import astropy.units as u
import matplotlib.patches as mpatches
import numpy as np
try:
    import pylab as plt
except ModuleNotFoundError:
    pass

eq_frames = ['icrs', 'fk5', 'fk4']

def sector_mask(shape, centre, radius, angle_range):
    """
    Return a boolean mask for a circular sector. The start/stop angles in  
    `angle_range` should be given in clockwise order.
    """

    x,y = np.ogrid[:shape[0],:shape[1]]
    cx,cy = centre
    tmin,tmax = np.deg2rad(angle_range)

    # ensure stop angle > start angle
    if tmax < tmin:
            tmax += 2*np.pi

    # convert cartesian --> polar coordinates
    r2 = (x-cx)*(x-cx) + (y-cy)*(y-cy)
    theta = np.arctan2(x-cx, y-cy) - tmin

    # wrap angles between 0 and 2*pi
    theta %= (2*np.pi)

    # circular mask
    circmask = r2 <= radius*radius

    # angular mask
    anglemask = theta <= (tmax-tmin)

    return circmask*anglemask

def ellipse_mask(shape, x0, y0, bmaj, bmin, angle):
    """
    """
    
    x,y = np.mgrid[:shape[0],:shape[1]]
    
    cos_angle = np.cos(np.radians(180. - angle))
    sin_angle = np.sin(np.radians(180. - angle))
    
    # Shift the indexes.
    xc = x - x0
    yc = y - y0
    
    # Rotate the indexes.
    xct = xc * cos_angle - yc * sin_angle
    yct = xc * sin_angle - yc * cos_angle
    
    mask = (xct/bmaj)**2. + (yct/bmin)**2. <= 1.
    
    return mask

def add_radius_units(value, units):
    """
    """
    
    if units == 'd':
        r = value*u.deg
    if units == 'm':
        r = value*u.arcmin
    if units == 's':<commit_msg>Fixing elliptical region extraction<commit_after>#!/usr/bin/env python3

"""
Extracts a spectrum from a region in a spectral cube.
The pixels will be averaged spatially inside the given region.
The region must be specified as shape,coords,frame,parameters
Region can be, point, box, circle or ellipse. Also, CASA region files are supported.
Coords can be pix or sky.
Frames supported include: fk4, fk5, icrs, and gal.
Parameters, for point the coordinates of the point, e.g. point,pix,,512,256.
For box the bottom left corner and top right corner coordinates, e.g. box,pix,,256,256,512,512.
For circle the center of the circle and the radius, e.g. circle,pix,,512,256,10.
For ellipse, its center, major and minor axes and the angle in degrees, e.g., 
ellipse,sky,fk5,13h05m27.462s,-49d28m06.547s,0.345293s,0.493276s,0
If the coordinates are given in sky values, then the units must included, e.g. 
point,sky,12h,2d will extract the spectrum for the pixel located at RA 12 hours and DEC 2 degrees.
For circle the radius units in sky coordinates can be either d for degrees, m for arcminutes or 
s for arcseconds. The conversion will use the scale in the RA direction of the cube, i.e. will 
use CDELT1 to convert from angular units to pixels.
"""

import os
import matplotlib as mpl
havedisplay = "DISPLAY" in os.environ
if not havedisplay:
    mpl.use('Agg')

import sys
import re
import argparse
import logging

from astropy.io import fits
from astropy.table import Table
from astropy.io import ascii
from astropy.wcs import WCS, WCSSUB_SPECTRAL
from astropy.coordinates import SkyCoord

import crrlpy.imtools as ci
import astropy.units as u
import matplotlib.patches as mpatches
import numpy as np
try:
    import pylab as plt
except ModuleNotFoundError:
    pass

eq_frames = ['icrs', 'fk5', 'fk4']

def sector_mask(shape, centre, radius, angle_range):
    """
    Return a boolean mask for a circular sector. The start/stop angles in  
    `angle_range` should be given in clockwise order.
    """

    x,y = np.ogrid[:shape[0],:shape[1]]
    cx,cy = centre
    tmin,tmax = np.deg2rad(angle_range)

    # ensure stop angle > start angle
    if tmax < tmin:
            tmax += 2*np.pi

    # convert cartesian --> polar coordinates
    r2 = (x-cx)*(x-cx) + (y-cy)*(y-cy)
    theta = np.arctan2(x-cx, y-cy) - tmin

    # wrap angles between 0 and 2*pi
    theta %= (2*np.pi)

    # circular mask
    circmask = r2 <= radius*radius

    # angular mask
    anglemask = theta <= (tmax-tmin)

    return circmask*anglemask

def ellipse_mask(shape, x0, y0, bmaj, bmin, angle):
    """
    """
    
    x,y = np.mgrid[:shape[0],:shape[1]]
    
    cos_angle = np.cos(np.radians(180. - angle))
    sin_angle = np.sin(np.radians(180. - angle))
    
    # Shift the indexes.
    xc = x - x0
    yc = y - y0
    
    # Rotate the indexes.
    xct = xc * cos_angle - yc * sin_angle
    yct = xc * sin_angle - yc * cos_angle
    
    mask = ( (xct/bmaj)**2. + (yct/bmin)**2. ) <= 1.
    
    return mask

def add_radius_units(value, units):
    """
    """
    
    if units == 'd':
        r = value*u.deg
    if units == 'm':
        r = value*u.arcmin
    if units == 's':
{'fraction_non_alphanumeric': 0.10307643514113543, 'fraction_numerical': 0.030764351411354265, 'mean_word_length': 3.0901426718547342, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '733737', 'n_tokens_mistral': 2323, 'n_tokens_neox': 2148, 'n_words': 770}