diff --git a/ds9_helpers/select_facets.py b/ds9_helpers/select_facets.py new file mode 100644 index 00000000..2458f505 --- /dev/null +++ b/ds9_helpers/select_facets.py @@ -0,0 +1,132 @@ +import pandas as pd +import numpy as np +from math import radians, cos, sin, sqrt, atan2 +from scipy.spatial import Voronoi, voronoi_plot_2d +import matplotlib.cm as cm +import matplotlib.pyplot as plt + + +def angular_distance(ra1_deg, dec1_deg, ra2_deg, dec2_deg): + """ + Function to calculate angular distance between two points with RA/DEC + + :param ra1_deg: RA first source + :param dec1_deg: DEC first source + :param ra2_deg: RA second source + :param dec2_deg: DEC second source + + :return: distance in degrees + + """ + # Convert degrees to radians + ra1 = radians(ra1_deg) + dec1 = radians(dec1_deg) + ra2 = radians(ra2_deg) + dec2 = radians(dec2_deg) + + # Haversine-like formula for angular distance in degrees + delta_ra = ra2 - ra1 + delta_dec = dec2 - dec1 + + a = sin(delta_dec / 2) ** 2 + cos(dec1) * cos(dec2) * sin(delta_ra / 2) ** 2 + c = 2 * atan2(sqrt(a), sqrt(1 - a)) + + # Convert back to degrees + distance_deg = np.degrees(c) + return distance_deg + + +def remove_close_pairs(csv_file, dist_offset=0.1): + """ + Remove close pairs in CSV and sort by scalarphasediff score + + Args: + csv_file: CSV file with sources RA/DEC and spd_score + dist_offset: min offset in degrees + + Returns: + cleaned df + """ + + # Load the CSV file + df = pd.read_csv(csv_file).sort_values('spd_score') + + # Assuming the CSV file has columns 'RA' and 'DEC' in degrees + ra = df['RA'].values + dec = df['DEC'].values + + # List to store pairs and their distances + pairs = [] + + # Calculate pairwise distances + for i in range(len(ra)): + for j in range(i + 1, len(ra)): + dist = angular_distance(ra[i], dec[i], ra[j], dec[j]) + if dist < dist_offset: + pairs.append((i, j, dist)) + + # Convert to a DataFrame to display + duplicate_df = pd.DataFrame(pairs, columns=['Index 1', 'Index 2', 'Distance (degrees)']) + + for idx in duplicate_df['Index 2'].values[::-1]: + df.drop(idx, inplace=True) + + return df + + +def plot_voronoi(csv_path, ra_column='RA', dec_column='DEC', score_column='spd_score'): + """ + Plots a Voronoi diagram where the points (RA/DEC) are colored based on the values in the score_column. + + Parameters: + csv_path (str): Path to the CSV file. + ra_column (str): Column name for RA values. + dec_column (str): Column name for DEC values. + score_column (str): Column name for the score values used for coloring the points. + """ + # Load the CSV data + data = pd.read_csv(csv_path).sort_values(score_column) + + # Extract RA, DEC, and spd_scores columns + ra = data[ra_column].values + dec = data[dec_column].values + spd_scores = data[score_column].values + + # Combine RA and DEC into coordinate pairs + points = np.column_stack((ra, dec)) + + # Compute Voronoi tessellation + vor = Voronoi(points) + + # Create a plot + fig, ax = plt.subplots(figsize=(10, 8)) + + # Plot the Voronoi diagram without filling the regions + voronoi_plot_2d(vor, ax=ax, show_vertices=False, line_colors='black', line_width=2, line_alpha=0.6) + + # Normalize spd_scores for coloring + norm = plt.Normalize(vmin=0, vmax=2) + cmap = cm.viridis + + # Scatter plot the points, colored by spd_scores + sc = ax.scatter(ra[1:], dec[1:], c=spd_scores[1:], cmap=cmap, norm=norm, edgecolor='black', s=200, zorder=2) + + # Highlight the first point with a red star + ax.scatter(ra[0], dec[0], color='red', marker='*', s=600, zorder=3, edgecolor='black') + + + ax.tick_params(labelsize=16) + + # Add a colorbar + # Add a colorbar and set font sizes + cbar = plt.colorbar(sc) + cbar.set_label('$\hat{\sigma}_{c}$ (rad)', fontsize=16) # Set the font size for the colorbar label + cbar.ax.tick_params(labelsize=16) # Set the font size for colorbar ticks + plt.gca().invert_xaxis() + + # Set labels and title + ax.set_xlabel('RA (degrees)', fontsize=16) + ax.set_ylabel('DEC (degrees)', fontsize=16) + + plt.tight_layout() + plt.savefig('voronoi_calibrators.png', dpi=150) \ No newline at end of file diff --git a/fits_helpers/cropy_2048.py b/fits_helpers/cropy_2048.py new file mode 100644 index 00000000..7e5c5ab1 --- /dev/null +++ b/fits_helpers/cropy_2048.py @@ -0,0 +1,57 @@ +from astropy.io import fits +import numpy as np + + +def crop_fits_image(input_filename, output_filename, center=None): + """ + Crops a FITS image to 2048x2048 pixels. + + Parameters: + - input_filename: str, path to the input FITS file. + - output_filename: str, path to the output cropped FITS file. + - center: tuple (x_center, y_center), the pixel coordinates of the center of the crop. + If None, the image will be cropped from the center of the input image. + """ + # Open the FITS file + with fits.open(input_filename) as hdul: + data = hdul[0].data.squeeze() + header = hdul[0].header + + # Get the dimensions of the image + y_size, x_size = data.shape + + # Determine the center of the crop + if center is None: + x_center = x_size // 2 + y_center = y_size // 2 + else: + x_center, y_center = center + + # Calculate the starting and ending indices for the crop + x_start = int(x_center - 1024) + x_end = int(x_center + 1024) + y_start = int(y_center - 1024) + y_end = int(y_center + 1024) + + # Ensure the indices are within the bounds of the image + x_start = max(0, x_start) + y_start = max(0, y_start) + x_end = min(x_size, x_end) + y_end = min(y_size, y_end) + + # Crop the image data + cropped_data = data[y_start:y_end, x_start:x_end] + + # Update header to reflect new image size + header['NAXIS1'] = cropped_data.shape[1] + header['NAXIS2'] = cropped_data.shape[0] + + # Adjust the reference pixel if present + if 'CRPIX1' in header: + header['CRPIX1'] -= x_start + if 'CRPIX2' in header: + header['CRPIX2'] -= y_start + + # Write the cropped image to the output file + hdu = fits.PrimaryHDU(data=np.array([[cropped_data]]), header=header) + hdu.writeto(output_filename, overwrite=True) diff --git a/ms_helpers/other/calc_uvw.py b/ms_helpers/other/calc_uvw.py index 543a61e5..6d0c7930 100644 --- a/ms_helpers/other/calc_uvw.py +++ b/ms_helpers/other/calc_uvw.py @@ -1,296 +1,98 @@ -from __future__ import print_function, division +from casacore.tables import table import numpy as np -from astropy.time import Time, TimeDelta -import matplotlib.pyplot as plt +from astropy.coordinates import EarthLocation, SkyCoord, GCRS +from astropy.time import Time import astropy.units as u -from astropy.coordinates import (SkyCoord, FK5, ITRS, CIRS, EarthLocation, CartesianRepresentation, - solar_system_ephemeris, get_body_barycentric_posvel) -from PyAstronomy import pyasl -import datetime +from astropy.utils import iers -def calculate_uvw(X1, Y1, Z1, X2, Y2, Z2, time_str, alpha, delta): - """ - Calculate UVW coordinates for a given observation. - - Parameters: - X1, Y1, Z1 : float : Geocentric coordinates of antenna 1 - X2, Y2, Z2 : float : Geocentric coordinates of antenna 2 - time_str : str : Observation time in UTC (ISO format string) - alpha : float : Right Ascension of the source in degrees - delta : float : Declination of the source in degrees - latitude : float : Latitude of the observation location in degrees - longitude : float : Longitude of the observation location in degrees - - Returns: - uvw : numpy.ndarray : UVW coordinates - """ - - def midpoint_to_lat_lon(X1, Y1, Z1, X2, Y2, Z2): - # Calculate midpoint coordinates - X_m = (X1 + X2) / 2 - Y_m = (Y1 + Y2) / 2 - Z_m = (Z1 + Z2) / 2 - - # Create a Cartesian representation of the midpoint - mid_cartesian = CartesianRepresentation(X_m * u.m, Y_m * u.m, Z_m * u.m) - - # Convert to EarthLocation to get latitude and longitude - mid_location = EarthLocation.from_geocentric(mid_cartesian.x, mid_cartesian.y, mid_cartesian.z) - - return mid_location.lat.deg, mid_location.lon.deg - - # Baseline vector - B = np.array([X2 - X1, Y2 - Y1, Z2 - Z1]) - - # Define observer location - latitude, longitude = midpoint_to_lat_lon(X1, Y1, Z1, X2, Y2, Z2) - location = EarthLocation(lat=latitude, lon=longitude) - - # Convert to LST - observation_time = Time(time_str, location=location) - lst = observation_time.sidereal_time('mean', longitude=longitude).deg - - # Calculate the hour angle - H = lst - alpha - if H < 0: - H += 360.0 - - # Convert degrees to radians - H_rad = np.radians(H) - delta_rad = np.radians(delta) - - # UVW Transformation matrix - transformation_matrix = np.array([ - [np.sin(H_rad), np.cos(H_rad), 0], - [-np.sin(delta_rad) * np.cos(H_rad), np.sin(delta_rad) * np.sin(H_rad), np.cos(delta_rad)], - [np.cos(delta_rad) * np.cos(H_rad), -np.cos(delta_rad) * np.sin(H_rad), np.sin(delta_rad)] - ]) - - # Compute UVW coordinates - uvw = np.dot(transformation_matrix, B) - return uvw - - -# def aberration_shift(timestamp, ra, dec): -# -# """ FROM: https://pyastronomy.readthedocs.io/en/latest/pyaslDoc/aslDoc/aberration.html """ -# -# # Convert calendar date to JD using the datetime package -# dt = datetime.datetime.strptime(timestamp, "%Y-%m-%dT%H:%M:%S") -# jd = pyasl.jdcnv(dt) -# -# # Specify RA and DEC -# ra = np.atleast_1d(ra) -# dec = np.atleast_1d(dec) -# -# # Get change in RA and DEC due to annual aberration -# res = pyasl.co_aberration(np.repeat(jd, ra.size), ra, dec) -# -# return res[0][0]*3600, res[1][0]*3600 - - -def aberration_correction(timestamp, ra, dec): - """ - Calculate the aberration corrected RA and DEC for given input RA, DEC and timestamp. - - Parameters: - ra (float): Right Ascension in degrees - dec (float): Declination in degrees - timestamp (str): Timestamp in ISO format (e.g., '2024-07-17T00:00:00') - - Returns: - corrected_ra (float): Aberration corrected Right Ascension in degrees - corrected_dec (float): Aberration corrected Declination in degrees - """ - # Convert timestamp to Julian Date - time = Time(timestamp, format='isot', scale='utc') - - # Initial coordinates - sky_coord = SkyCoord(ra=ra, dec=dec, unit='deg') +# Ensure the latest Earth Orientation Parameters are used +iers.conf.auto_download = True - # Calculate the position and velocity of Earth - with solar_system_ephemeris.set('builtin'): - pos, vel = get_body_barycentric_posvel('earth', time) - - # Speed of light in AU/day - c = 173.144632674240 - - # Aberration correction - ra_rad = np.deg2rad(ra) - dec_rad = np.deg2rad(dec) - - # Earth's velocity components - vx = vel.x.value - vy = vel.y.value - vz = vel.z.value - - # Proper motion in RA and DEC - delta_ra = (-vx * np.sin(ra_rad) + vy * np.cos(ra_rad)) / c - delta_dec = (-vx * np.cos(ra_rad) * np.sin(dec_rad) - - vy * np.sin(ra_rad) * np.sin(dec_rad) - + vz * np.cos(dec_rad)) / c - - - return np.rad2deg(delta_ra) * 3600, np.rad2deg(delta_dec) * 3600 - - -def precession_nutation_shift(timestamp, ra, dec): - """ - Calculate the precession and nutation corrected RA and DEC for given input RA, DEC, and timestamp. - - Parameters: - ra (float): Right Ascension in degrees - dec (float): Declination in degrees - timestamp (str): Timestamp in ISO format (e.g., '2024-07-17T00:00:00') - - Returns: - corrected_ra (float): Precession and nutation corrected Right Ascension in degrees - corrected_dec (float): Precession and nutation corrected Declination in degrees - """ - # Convert timestamp to astropy Time object - t = Time(timestamp) - - # Create a SkyCoord object for the input RA and DEC - coord = SkyCoord(ra=ra * u.degree, dec=dec * u.degree, frame='icrs') - - # Convert the coordinates to the FK5 frame at the given epoch to account for precession - coord_fk5 = coord.transform_to(FK5(equinox=t)) - - # Convert the coordinates to the CIRS frame to account for nutation - coord_cirs = coord_fk5.transform_to(CIRS(obstime=t)) - - # Extract the corrected RA and DEC - corrected_ra = coord_cirs.ra.degree - corrected_dec = coord_cirs.dec.degree - - return (ra-corrected_ra)*3600, (dec-corrected_dec)*3600 - - -def calculate_uvw_shift(ra_deg, dec_deg, baseline_length, ra_shift_arcsec, dec_shift_arcsec): +def recalculate_uvw(ms_path): """ - Calculate the shift in UVW coordinates based on shifts in RA and DEC. + Recalculate UVW coordinates for a measurement set. Parameters: - ra_deg (float): Original Right Ascension in degrees - dec_deg (float): Original Declination in degrees - baseline_length (float): Baseline length in meters - ra_shift_arcsec (float): Shift in Right Ascension in arcseconds - dec_shift_arcsec (float): Shift in Declination in arcseconds + - ms_path (str): Path to the measurement set. Returns: - delta_u (float): Shift in U coordinate in meters - delta_v (float): Shift in V coordinate in meters - delta_w (float): Shift in W coordinate in meters + - None: The function updates the UVW coordinates in the MS in-place. """ - def ra_dec_to_radians(ra_shift_arcsec, dec_shift_arcsec): - """Convert RA and DEC shifts from arcseconds to radians.""" - ra_shift_rad = np.deg2rad(ra_shift_arcsec / 3600.0) - dec_shift_rad = np.deg2rad(dec_shift_arcsec / 3600.0) - return ra_shift_rad, dec_shift_rad - - def direction_cosines(ra_rad, dec_rad): - """Calculate direction cosines from RA and DEC in radians.""" - l = np.cos(dec_rad) * np.cos(ra_rad) - m = np.cos(dec_rad) * np.sin(ra_rad) - n = np.sin(dec_rad) - return l, m, n - - def uvw_shift(ra_rad, dec_rad, ra_shift_rad, dec_shift_rad, baseline_length): - """Calculate the shift in UVW coordinates based on RA and DEC shifts.""" - # Original direction cosines - l0, m0, n0 = direction_cosines(ra_rad, dec_rad) - - # Shifted RA and DEC - ra_shifted_rad = ra_rad + ra_shift_rad - dec_shifted_rad = dec_rad + dec_shift_rad - - # Shifted direction cosines - l1, m1, n1 = direction_cosines(ra_shifted_rad, dec_shifted_rad) - - # Difference in direction cosines - delta_l = l1 - l0 - delta_m = m1 - m0 - delta_n = n1 - n0 - - # Convert direction cosine shifts to UVW coordinate shifts - delta_u = baseline_length * delta_l - delta_v = baseline_length * delta_m - delta_w = baseline_length * delta_n - - return delta_u, delta_v, delta_w + # Open the measurement set in writable mode + with table(ms_path, readonly=False, ack=False) as ms: + # Read necessary columns from the MAIN table + time_col = ms.getcol('TIME') # Observation times in seconds since MJD epoch + antenna1_col = ms.getcol('ANTENNA1') # Indices of first antennas + antenna2_col = ms.getcol('ANTENNA2') # Indices of second antennas - # Convert shifts to radians - ra_shift_rad, dec_shift_rad = ra_dec_to_radians(ra_shift_arcsec, dec_shift_arcsec) + # Get antenna positions from the ANTENNA table + with table(ms_path + '::ANTENNA', ack=False) as ant_table: + ant_positions = ant_table.getcol('POSITION') # Shape: (n_antennas, 3) + n_antennas = ant_table.nrows() - # Convert RA and DEC to radians - ra_rad = np.deg2rad(ra_deg) - dec_rad = np.deg2rad(dec_deg) + # Get source direction from the FIELD table + with table(ms_path + '::FIELD', ack=False) as field_table: + phase_dir = field_table.getcol('PHASE_DIR') # Shape: (nfields, 1, 2) + ra = phase_dir[0, 0, 0] # Right Ascension in radians + dec = phase_dir[0, 0, 1] # Declination in radians - # Calculate UVW coordinate shifts - delta_u, delta_v, delta_w = uvw_shift(ra_rad, dec_rad, ra_shift_rad, dec_shift_rad, baseline_length) + # Convert times to astropy Time objects + times = Time(time_col / 86400.0, format='mjd', scale='utc') - return delta_u, delta_v, delta_w + # Convert antenna positions to EarthLocation objects + ant_locations = EarthLocation(x=ant_positions[:, 0] * u.m, + y=ant_positions[:, 1] * u.m, + z=ant_positions[:, 2] * u.m) + # Initialize arrays to hold GCRS positions + ant_gcrs = [] -# Example usage: -X1, Y1, Z1 = 3826896.235000, 460979.455000, 5064658.203000 -X2, Y2, Z2 = 3370271.657000, 712125.881000, 5349991.165000 -ra = 241.88 # Right Ascension in degrees -dec = 54.33 # Declination in degrees + # Transform antenna positions from ITRF to GCRS at each time + # Since the antennas are fixed on Earth, their positions in GCRS change over time due to Earth's rotation. -# Length of a sidereal day -sidereal_day = TimeDelta(23.9344696 * 3600, format='sec') # 23 hours, 56 minutes, 4 seconds + print("Transforming antenna positions to GCRS frame...") + for ant_loc in ant_locations: + itrs = ant_loc.get_itrs(obstime=times[0]) + gcrs = itrs.transform_to(GCRS(obstime=times[0])) + ant_gcrs.append(gcrs.cartesian) + ant_gcrs = np.array([[gcrs.x.value, gcrs.y.value, gcrs.z.value] for gcrs in ant_gcrs]) # Shape: (n_antennas, 3) -for n, timez in enumerate(["2018-11-26T07:13:43", "2020-05-24T19:20:26", "2020-11-14T08:11:00", "2021-05-13T19:41:00"]): + # Compute baselines in GCRS frame + bl_gcrs = ant_gcrs[antenna2_col] - ant_gcrs[antenna1_col] # Shape: (nrows, 3) - dra_aber, ddec_aber = aberration_correction(timez, ra, dec) - dra_prec, ddec_prec = precession_nutation_shift(timez, ra, dec) - # if n>=1: - # ddec_aber-=40 - du, dv, dw = calculate_uvw_shift(ra, dec, 600_000, dra_prec + dra_aber, ddec_prec + ddec_aber) - print(timez) - print(dra_prec, ddec_prec) - print(dra_aber, ddec_aber) + # Define the source direction in GCRS frame + source_icrs = SkyCoord(ra=ra * u.rad, dec=dec * u.rad, frame='icrs') + source_gcrs = source_icrs.transform_to(GCRS(obstime=times[0])) # Use the first time for transformation + # Get the unit vector pointing to the source in GCRS frame + src_unit_vector = np.array([source_gcrs.cartesian.x.value, + source_gcrs.cartesian.y.value, + source_gcrs.cartesian.z.value]) # Shape: (3,) - if n==0: - ref_dra_prec, ref_ddec_prec = dra_prec, ddec_prec - ref_dra_aber, ref_ddec_aber = dra_aber, ddec_aber - ref_du, ref_dv, ref_dw = du, dv, dw + # Create orthonormal basis (u_vec, v_vec, w_vec) where w_vec is towards the source + w_vec = src_unit_vector / np.linalg.norm(src_unit_vector) # Normalize w_vec - # else: - # print(dra_prec-ref_dra_prec, ddec_prec-ref_ddec_prec) - # print(dra_aber-ref_dra_aber, ddec_aber-ref_ddec_aber) - # print((dra_prec-ref_dra_prec)-(dra_aber-ref_dra_aber), (ddec_prec-ref_ddec_prec)-(ddec_aber-ref_ddec_aber)) + # Choose an arbitrary vector not parallel to w_vec to construct u_vec + arbitrary_vector = np.array([0, 0, 1]) + if np.allclose(w_vec, arbitrary_vector): + arbitrary_vector = np.array([0, 1, 0]) + # Compute u_vec = cross(arbitrary_vector, w_vec) + u_vec = np.cross(arbitrary_vector, w_vec) + u_vec /= np.linalg.norm(u_vec) # Normalize u_vec + # Compute v_vec = cross(w_vec, u_vec) + v_vec = np.cross(w_vec, u_vec) - print(du-ref_du, dv-ref_dv) - if n==0: - plt.scatter([(du-ref_du)/1000+518.45-i*0.03 for i in range(10)], - [(dv-ref_dv)/1000+212.55 + 0.3*i for i in range(10)], label=timez, marker="*") - else: - plt.scatter([(du-ref_du)/1000-0.033+518.45-i*0.03 for i in range(10)], - [(dv-ref_dv)/1000-0.033+212.55 + 0.3*i for i in range(10)], label=timez, marker="*") - print() + # Now project the baselines onto the (u_vec, v_vec, w_vec) basis + uvw = np.zeros_like(bl_gcrs) + uvw[:, 0] = np.dot(bl_gcrs, u_vec) + uvw[:, 1] = np.dot(bl_gcrs, v_vec) + uvw[:, 2] = np.dot(bl_gcrs, w_vec) - # # pos = np.array( - # # [calculate_uvw(X1, Y1, Z1, X2, Y2, Z2, - # # (Time(timez) + i*4/(23.9344696 * 3600)).iso, - # # ra-(dra_prec+dra_aber)/3600, - # # dec-(ddec_prec+ddec_aber)/3600) - # # for i in range(int(3600))]) - # pos = np.array( - # [calculate_uvw(X1, Y1, Z1, X2, Y2, Z2, - # (Time(timez) + i*8/(23.9344696 * 3600)).iso, - # ra+(dra_prec+dra_aber)/3600, - # dec+(ddec_prec+ddec_aber)/3600) - # for i in range(int(3600))]) - # plt.scatter(pos[:, 0]/1000, pos[:, 1]/1000, label=timez, s=30, marker="*") + # Update the UVW column in the measurement set + ms.putcol('UVW', uvw) -plt.xlim(518.1, 518.6) -plt.ylim(212.4, 215.1) -plt.legend() -plt.show() + print("UVW recalculation completed.") diff --git a/ms_helpers/other/doppler_correction.py b/ms_helpers/other/doppler_correction.py deleted file mode 100644 index 33afe936..00000000 --- a/ms_helpers/other/doppler_correction.py +++ /dev/null @@ -1,101 +0,0 @@ -import numpy as np -import casacore.tables as tables -import casacore.measures as measures -import casacore.quanta as quanta - -def calculate_radial_velocity(observation_time, ant_position, source_coords, frame='LSRK'): - me = measures.measures() - qa = quanta.quanta() - - # Convert observation time to casacore quantity - time_quantity = qa.quantity(observation_time, 's') - - # Set up frames - me.doframe(me.epoch('UTC', time_quantity)) - me.doframe(me.position('ITRF', *[qa.quantity(pos, 'm') for pos in ant_position])) - me.doframe(me.direction('J2000', *[qa.quantity(coord, 'rad') for coord in source_coords])) - - # Calculate radial velocity in the specified frame - doppler = me.doppler('RADIO', frame) - radial_velocity = doppler['m0']['value'] # in m/s - return radial_velocity - -def correct_doppler_shift_casacore(msname, restfreq, outms, frame='LSRK'): - """ - Corrects Doppler shifts in a Measurement Set using casacore. - - Parameters: - msname (str): Name of the input Measurement Set. - restfreq (str): Rest frequency of the observed line (e.g., '1420.40575177MHz'). - outms (str): Name of the output Measurement Set with Doppler correction applied. - frame (str): Reference frame for the Doppler correction (default is 'LSRK'). - - Returns: - None - """ - # Initialize casacore tools - me = measures.measures() - qa = quanta.quanta() - - # Open the Measurement Set - ms = tables.table(msname, readonly=False) - - # Get the observation time and antenna positions - obs_table = tables.table(msname + '/OBSERVATION') - obs_time = obs_table.getcol('TIME_RANGE')[0][0] - obs_table.close() - - ant_table = tables.table(msname + '/ANTENNA') - ant_positions = ant_table.getcol('POSITION') - ant_table.close() - - # Get the field direction - field_table = tables.table(msname + '/FIELD') - phase_dir = field_table.getcol('PHASE_DIR')[0] - field_table.close() - - # Convert the rest frequency to a casacore quantity - restfreq_qa = qa.quantity(restfreq) - - # Prepare to calculate Doppler shifts - me.doframe(me.epoch('UTC', qa.quantity(obs_time, 's'))) - avg_doppler = [] - - for i in range(len(ant_positions)): - ant_pos = me.position('ITRF', - qa.quantity(ant_positions[i][0], 'm'), - qa.quantity(ant_positions[i][1], 'm'), - qa.quantity(ant_positions[i][2], 'm')) - me.doframe(ant_pos) - me.doframe(me.direction('J2000', - qa.quantity(phase_dir[0], 'rad'), - qa.quantity(phase_dir[1], 'rad'))) - doppler = me.measure(me.doppler('RADIO', frame, restfreq_qa), frame) - avg_doppler.append(doppler['m0']['value']) - - avg_doppler_shift = np.mean(avg_doppler) - - # Apply the Doppler correction factor - doppler_factor = 1 + avg_doppler_shift / 299792.458 # speed of light in km/s - - # Update the frequencies in the SPW table - spw_table = tables.table(msname + '/SPECTRAL_WINDOW', readonly=False) - chan_freqs = spw_table.getcol('CHAN_FREQ') - - corrected_chan_freqs = chan_freqs / doppler_factor - spw_table.putcol('CHAN_FREQ', corrected_chan_freqs) - spw_table.close() - - # Copy the Measurement Set to the output file - ms.copy(outms) - ms.close() - - print(f"Doppler correction applied. Output MS saved as {outms}") - - -if __name__ == '__main__': - # Example usage: - msname = 'your_measurement_set.ms' - restfreq = '1420.40575177MHz' - outms = 'corrected_measurement_set.ms' - correct_doppler_shift_casacore(msname, restfreq, outms) \ No newline at end of file diff --git a/ms_helpers/other/doppler_shift.py b/ms_helpers/other/doppler_shift.py index a3e6108c..dae0c3a1 100644 --- a/ms_helpers/other/doppler_shift.py +++ b/ms_helpers/other/doppler_shift.py @@ -1,38 +1,44 @@ -from astropy.coordinates import SkyCoord, EarthLocation, AltAz +from casacore.tables import table +from astropy.coordinates import SkyCoord, EarthLocation from astropy.time import Time import astropy.units as u from scipy.constants import c -# Define the observing frequency -observing_frequency = 140 * u.MHz +def calculate_doppler_shifts(ms_path, instrument="LOFAR"): + """ + Calculate Doppler shifts for an MS. -# Define the location of the observatory (example: Very Large Array) -observatory_location = EarthLocation.of_site('LOFAR') + Parameters: + - ms_path (str): Path to the MS file. + - instrument (str): Instrument name. -# Define the source coordinates -source = SkyCoord(ra=-2.04639855*u.rad, dec=0.95905842*u.rad, frame='icrs') + Returns: + - result: Doppler shift in kHz + """ -for t in ['2018-11-26', '2020-05-24', '2020-11-14', '2021-5-13']: - print(t) + # Get observation time + with table(ms_path+"::POINTING", ack=False) as ms: + obs_time_seconds = ms.getcol('TIME_ORIGIN')[0] # Extract observation times in seconds + ra, dec = ms.getcol("DIRECTION")[0][0] - # Define the observation times - time = Time(t) + # Convert seconds to MJD + obs_time_mjd = obs_time_seconds / 86400.0 + time = Time(obs_time_mjd, format='mjd', scale='utc') - # Calculate the AltAz frame for the observation times - altaz = AltAz(location=observatory_location, obstime=time) + # Get ref frequency + with table(ms_path+"::SPECTRAL_WINDOW", ack=False) as ms: + observing_frequency = ms.getcol("REF_FREQUENCY")[0] * u.Hz - # Transform the source coordinates to the AltAz frame - altaz_coord = source.transform_to(altaz) + # Define the location of the observatory + observatory_location = EarthLocation.of_site(instrument) - # Calculate the radial velocity due to Earth's motion - rv = altaz_coord.radial_velocity_correction() - - # Calculate the Doppler shift for each observation - doppler_shift = observing_frequency * (rv / c / u.m*u.s) + # Define the source coordinates + source = SkyCoord(ra=ra * u.rad, dec=dec * u.rad, frame='icrs') - print(f"Doppler shift for observation 1: {doppler_shift.to(u.kHz)}") + # Calculate the radial velocity due to Earth's motion + rv = source.radial_velocity_correction(obstime=time, location=observatory_location).to(u.m / u.s) - # Adjust the observed frequencies before stacking - adjusted_frequency = observing_frequency + doppler_shift + # Calculate the Doppler shift for the observation + doppler_shift = observing_frequency * (rv / c) - print(f"Adjusted frequency for observation 1: {adjusted_frequency.to(u.MHz)}") + return doppler_shift.to(u.kHz)