From fade757c0c63570f50d6d821347c520861bea6cc Mon Sep 17 00:00:00 2001 From: HrishikeshShetgaonkar Date: Tue, 25 Aug 2026 09:10:21 +0200 Subject: [PATCH 1/3] added linear polarization plotting --- src/plot_field/plot_field.py | 201 ++++++++++++++++++++++++++++++++++- 1 file changed, 200 insertions(+), 1 deletion(-) diff --git a/src/plot_field/plot_field.py b/src/plot_field/plot_field.py index b1073d3..3bf53af 100644 --- a/src/plot_field/plot_field.py +++ b/src/plot_field/plot_field.py @@ -588,6 +588,7 @@ def make_plot( lbcs_catalogue, targRA=None, targDEC=None, + polarized_sources=None, nchan=16, av_time=1, outdir=".", @@ -684,6 +685,7 @@ def make_plot( label="LBCS Sources", ) + c = SphericalCircle( centre_coord, 1.5 * u.deg, @@ -693,6 +695,22 @@ def make_plot( ) ax.add_patch(c) + # plot polarized sources + polint = np.asarray(polarized_sources["polint"], dtype=float) + marker_size = 20 + 100 * polint / np.nanmax(polint) #20 is just a factor to get bigger markers + ax.scatter( + polarized_sources["RA_pol"], + polarized_sources["Dec_pol"], + transform=ax.get_transform("fk5"), + marker="D", + s=marker_size, + facecolors="none", + edgecolors="magenta", + linewidths=1.5, + label="Polarized Sources", + zorder=5, + ) + ### Calculate and plot smallest distance to target dist_ids, dist = smallest_distance(RA, DEC, lbcs) @@ -895,6 +913,138 @@ def ps_match(file): t.write(file, overwrite=True) +def _parse_angular_radius(radius): + """Return an angular radius as an Astropy Quantity. + """ + if isinstance(radius, u.Quantity): + return radius.to(u.deg) + + if isinstance(radius, (int, float, np.integer, np.floating)): + return float(radius) * u.deg + + radius = str(radius).strip() + try: + quantity = u.Quantity(radius) + if quantity.unit == u.dimensionless_unscaled: + quantity = float(quantity.value) * u.deg + return quantity.to(u.deg) + except (ValueError, TypeError, u.UnitConversionError) as exc: + raise ValueError( + f"Could not parse RM-grid radius '{radius}'. " + "Use values such as '1.2deg' or '30arcmin'." + ) from exc + + +def rmgrid_match( + delay_cals_file, + rmgrid_catalogue, + field_ra, + field_dec, + field_radius="1.2deg", + pol_radius="0.3deg" +): + """find polarized sources within given RA,DEC + """ + if not os.path.isfile(rmgrid_catalogue): + raise FileNotFoundError(f"RM-grid catalogue not found: {rmgrid_catalogue}") + + rmgrid = Table.read(rmgrid_catalogue) + + rm_ra = "RA_pol" + rm_dec = "Dec_pol" + + pol_columns = ["RA_pol", "Dec_pol", "RM", "RM_err", "polint", "I_144MHz", "fracpol"] + missing = [name for name in pol_columns if name not in rmgrid.colnames] + if missing: + raise KeyError("Missing required RM-grid column(s): " + ", ".join(missing)) + + radius = _parse_angular_radius(field_radius) + field_centre = SkyCoord( + float(field_ra), float(field_dec), unit="deg", frame="icrs" + ) + rm_coords_all = SkyCoord( + rmgrid[rm_ra], rmgrid[rm_dec], unit="deg", frame="icrs" + ) + in_field = rm_coords_all.separation(field_centre) <= radius + polarized_sources = rmgrid[in_field] + + print( + f"Found {len(polarized_sources)} polarized sources within " + f"{radius.to_value(u.deg):.4g} deg " + f"({radius.to_value(u.arcmin):.4g} arcmin) of " + f"RA={float(field_ra):.6f}, DEC={float(field_dec):.6f}" + ) + + # Now write the data for nearest ones to delay cal. + delay_cals = Table.read(delay_cals_file) + search_radius = _parse_angular_radius(pol_radius) + + # Initialize output arrays + n_pol_sources = np.zeros(len(delay_cals), dtype=int) + nearest_pol_dist = np.full(len(delay_cals), np.nan) + nearest_polint = np.full(len(delay_cals), np.nan) + nearest_fracpol = np.full(len(delay_cals), np.nan) + nearest_pol_ra = np.full(len(delay_cals), np.nan) + nearest_pol_dec = np.full(len(delay_cals), np.nan) + + # find polarized sources near the delay cals + if len(polarized_sources) > 0: + pol_coords = SkyCoord( + polarized_sources[rm_ra], + polarized_sources[rm_dec], + unit="deg", + frame="icrs", + ) + + for i, delay_cal in enumerate(delay_cals): + + cal_coord = SkyCoord( + float(delay_cal["RA"]), + float(delay_cal["DEC"]), + unit="deg", + frame="icrs", + ) + + # Distances from this delay calibrator to all polarized sources + separations = cal_coord.separation(pol_coords) + + within_radius = separations <= search_radius + n_pol_sources[i] = np.count_nonzero(within_radius) + + if n_pol_sources[i] > 0: + idx_within = np.where(within_radius)[0] + nearest_idx = idx_within[np.argmin(separations[idx_within])] + + nearest_pol_dist[i] = separations[nearest_idx].deg + nearest_polint[i] = polarized_sources["polint"][nearest_idx] + nearest_fracpol[i] = polarized_sources["fracpol"][nearest_idx] + nearest_pol_ra[i] = polarized_sources[rm_ra][nearest_idx] + nearest_pol_dec[i] = polarized_sources[rm_dec][nearest_idx] + + # Add columns in delay calibrator table + new_columns = { + "n_pol": n_pol_sources, + "nearest_pol_dist": nearest_pol_dist, + "nearest_polint": nearest_polint, + "nearest_fracpol": nearest_fracpol, + "nearest_pol_ra": nearest_pol_ra, + "nearest_pol_dec": nearest_pol_dec, + } + + for name, values in new_columns.items(): + if name in delay_cals.colnames: + delay_cals[name] = values + else: + delay_cals.add_column(Column(values, name=name)) + + delay_cals.write( + delay_cals_file, + format="csv", + overwrite=True, + ) + + return polarized_sources + def gaia_quasar_match(file): from astroquery.vizier import Vizier @@ -1160,6 +1310,8 @@ def generate_catalogues( bright_limit_Jy=5.0, lotss_catalogue="lotss_catalogue.csv", lbcs_catalogue="lbcs_catalogue.csv", + rmgrid_catalogue="../../DR3_RMGrid_v2.fits", + pol2delay_radius='0.3deg', lotss_result_file="image_catalogue.csv", delay_cals_file="delay_calibrators.csv", match_tolerance=5.0, @@ -1398,6 +1550,14 @@ def generate_catalogues( gaia_quasar_match(delay_cals_file) ps_match(delay_cals_file) + polarized_sources = rmgrid_match( + delay_cals_file, + rmgrid_catalogue, + RATar, + DECTar, + im_radius, + pol_radius=pol2delay_radius + ) print("Assumed averaging - nchannels: %s; time averaging: %s" % (nchan, av_time)) make_plot( @@ -1409,6 +1569,7 @@ def generate_catalogues( targRA, targDEC, nchan=nchan, + polarized_sources=polarized_sources, av_time=av_time, outdir=outdir, ) @@ -1495,6 +1656,23 @@ def parse_args(): help="input file for LoTSS catalogue [will be downloaded if does not exist]", default="lotss_catalogue.csv", ) + parser.add_argument( + "--rmgrid_catalogue", + dest="rmgrid_catalogue", + type=str, + help="input file for rmgrid_catalog [will be downloaded if does not exist]", + default="../DR3_RMGrid_v2.fits", # ToDo: Set a default file path/url + ) + parser.add_argument( + "--pol2delay_radius", + dest="pol2delay_radius", + type=str, + help=( + "Radius around delay calibrator to choose polarized sources. " + "Examples: 0.3deg, 30arcmin [default: 0.3deg]" + ), + default="0.3deg", + ) parser.add_argument( "--lbcs_catalogue", dest="lbcs_catalogue", @@ -1620,6 +1798,8 @@ def main(): bright_limit_Jy=args.bright_limit_Jy, lotss_catalogue=args.lotss_catalogue, lbcs_catalogue=args.lbcs_catalogue, + rmgrid_catalogue=args.rmgrid_catalogue, + pol2delay_radius=args.pol2delay_radius, lotss_result_file=args.lotss_result_file, delay_cals_file=args.delay_cals_file, match_tolerance=args.match_tolerance, @@ -1672,6 +1852,23 @@ def main(): help="input file for LoTSS catalogue [will be downloaded if does not exist]", default="lotss_catalogue.csv", ) + parser.add_argument( + "--rmgrid_catalogue", + dest="rmgrid_catalogue", + type=str, + help="input file for rmgrid_catalog [will be downloaded if does not exist]", + default="../DR3_RMGrid_v2.fits", # ToDo: Set a default file path/url + ) + parser.add_argument( + "--pol2delay_radius", + dest="pol2delay_radius", + type=str, + help=( + "Radius around delay calibrators to select polarized sources. " + "Examples: 0.3deg, 30arcmin [default: 0.3deg]" + ), + default="0.3deg", + ) parser.add_argument( "--lbcs_catalogue", dest="lbcs_catalogue", @@ -1792,6 +1989,8 @@ def main(): lotss_catalogue=args.lotss_catalogue, lbcs_catalogue=args.lbcs_catalogue, lotss_result_file=args.lotss_result_file, + rmgrid_catalogue=args.rmgrid_catalogue, + pol2delay_radius=args.pol2delay_radius, delay_cals_file=args.delay_cals_file, match_tolerance=args.match_tolerance, image_limit_Jy=args.image_limit_Jy, @@ -1809,4 +2008,4 @@ def main(): ### TO DO LIST -# Add colour for quality of LBCS sources - Use fit or the compactness codes - Maybe not +# Add colour for quality of LBCS sources - Use fit or the compactness codes - Maybe not \ No newline at end of file From 01ff600854b8f30f42511d584566cc438355905e Mon Sep 17 00:00:00 2001 From: HrishikeshShetgaonkar Date: Thu, 27 Aug 2026 12:02:35 +0200 Subject: [PATCH 2/3] added pulsar plotting --- src/plot_field/plot_field.py | 143 +++++++++++++++++++++++++++++------ 1 file changed, 119 insertions(+), 24 deletions(-) diff --git a/src/plot_field/plot_field.py b/src/plot_field/plot_field.py index 17ef70f..129f2ea 100644 --- a/src/plot_field/plot_field.py +++ b/src/plot_field/plot_field.py @@ -5,7 +5,7 @@ import io import numpy as np import pyvo as vo -from astropy.table import Table, Column, hstack, unique +from astropy.table import Table, Column, hstack, unique, vstack import argparse # from lofarpipe.support.data_map import DataMap @@ -743,18 +743,36 @@ def make_plot( # plot polarized sources polint = np.asarray(polarized_sources["polint"], dtype=float) marker_size = 20 + 100 * polint / np.nanmax(polint) #20 is just a factor to get bigger markers + is_pulsar = polarized_sources["pol_type"] == "Pulsar" + is_rmgrid = ~is_pulsar + + # Normal polarized sources ax.scatter( - polarized_sources["RA_pol"], - polarized_sources["Dec_pol"], - transform=ax.get_transform("fk5"), - marker="D", - s=marker_size, - facecolors="none", - edgecolors="magenta", - linewidths=1.5, - label="Polarized Sources", - zorder=5, - ) + polarized_sources["RA_pol"][is_rmgrid], + polarized_sources["Dec_pol"][is_rmgrid], + transform=ax.get_transform("fk5"), + marker="D", + s=marker_size[is_rmgrid], + facecolors="none", + edgecolors="magenta", + linewidths=1.5, + label="Polarized Source", + zorder=5, + ) + + # Pulsars + ax.scatter( + polarized_sources["RA_pol"][is_pulsar], + polarized_sources["Dec_pol"][is_pulsar], + transform=ax.get_transform("fk5"), + marker="D", + s=marker_size[is_pulsar], + facecolors="none", + edgecolors="red", + linewidths=1.5, + label="Pulsars", + zorder=5, + ) ### Calculate and plot smallest distance to target dist_ids, dist = smallest_distance(RA, DEC, lbcs) @@ -985,12 +1003,14 @@ def _parse_angular_radius(radius): def rmgrid_match( delay_cals_file, rmgrid_catalogue, + pulsar_catalogue, field_ra, field_dec, field_radius="1.2deg", pol_radius="0.3deg" ): """find polarized sources within given RA,DEC + first for sources from DR3_RMGrid_v2.fits """ if not os.path.isfile(rmgrid_catalogue): raise FileNotFoundError(f"RM-grid catalogue not found: {rmgrid_catalogue}") @@ -1016,12 +1036,69 @@ def rmgrid_match( polarized_sources = rmgrid[in_field] print( - f"Found {len(polarized_sources)} polarized sources within " - f"{radius.to_value(u.deg):.4g} deg " - f"({radius.to_value(u.arcmin):.4g} arcmin) of " - f"RA={float(field_ra):.6f}, DEC={float(field_dec):.6f}" + f"Found {len(polarized_sources)} LoTSS polarized " + f"{radius.to_value(u.deg):.4g} deg " + f"({radius.to_value(u.arcmin):.4g} arcmin) of " + f"RA={float(field_ra):.6f}, DEC={float(field_dec):.6f}" ) + """ + now for pulsar sources from DR3_pulsars_known_v2.fits + """ + if not os.path.isfile(pulsar_catalogue): + raise FileNotFoundError(f"pulsar catalogue not found: {pulsar_catalogue}") + + pulsars = Table.read(pulsar_catalogue) + + pulsar_ra = "ra" + pulsar_dec = "dec" + + pol_columns = ["ra", "dec", "RM", "RM_err", "polint", "I_144MHz", "fracpol"] + missing = [name for name in pol_columns if name not in pulsars.colnames] + if missing: + raise KeyError("Missing required pulsar column(s): " + ", ".join(missing)) + + radius = _parse_angular_radius(field_radius) + field_centre = SkyCoord( + float(field_ra), float(field_dec), unit="deg", frame="icrs" + ) + pulsar_coords_all = SkyCoord( + pulsars[pulsar_ra], pulsars[pulsar_dec], unit="deg", frame="icrs" + ) + in_field = pulsar_coords_all.separation(field_centre) <= radius + in_field_pulsars = pulsars[in_field] + + print( + f"Found {len(in_field_pulsars)} pulsars " + f"{radius.to_value(u.deg):.4g} deg " + f"({radius.to_value(u.arcmin):.4g} arcmin) of " + f"RA={float(field_ra):.6f}, DEC={float(field_dec):.6f}" + ) + + # combine both catalogs + rmgrid_pol = polarized_sources[["RA_pol", "Dec_pol", "RM", "RM_err", "polint", "I_144MHz", "fracpol"]].copy() + rmgrid_pol["pol_type"] = ["RMGrid"] * len(rmgrid_pol) + # Pulsars + pulsar_pol = in_field_pulsars[ + ["ra", "dec", "RM", "RM_err", "polint", "I_144MHz", "fracpol"] + ].copy() + + pulsar_pol.rename_column("ra", "RA_pol") + pulsar_pol.rename_column("dec", "Dec_pol") + pulsar_pol["pol_type"] = ["Pulsar"] * len(pulsar_pol) + + # Combine them + if len(rmgrid_pol) > 0 and len(pulsar_pol) > 0: + all_polarized_sources = vstack( + [rmgrid_pol, pulsar_pol], + join_type="exact", + metadata_conflicts="silent", + ) + elif len(rmgrid_pol) <= 0: + all_polarized_sources = pulsar_pol + elif len(pulsar_pol) <= 0: + all_polarized_sources = rmgrid_pol + # Now write the data for nearest ones to delay cal. delay_cals = Table.read(delay_cals_file) search_radius = _parse_angular_radius(pol_radius) @@ -1035,10 +1112,10 @@ def rmgrid_match( nearest_pol_dec = np.full(len(delay_cals), np.nan) # find polarized sources near the delay cals - if len(polarized_sources) > 0: + if len(all_polarized_sources) > 0: pol_coords = SkyCoord( - polarized_sources[rm_ra], - polarized_sources[rm_dec], + all_polarized_sources[rm_ra], + all_polarized_sources[rm_dec], unit="deg", frame="icrs", ) @@ -1063,10 +1140,10 @@ def rmgrid_match( nearest_idx = idx_within[np.argmin(separations[idx_within])] nearest_pol_dist[i] = separations[nearest_idx].deg - nearest_polint[i] = polarized_sources["polint"][nearest_idx] - nearest_fracpol[i] = polarized_sources["fracpol"][nearest_idx] - nearest_pol_ra[i] = polarized_sources[rm_ra][nearest_idx] - nearest_pol_dec[i] = polarized_sources[rm_dec][nearest_idx] + nearest_polint[i] = all_polarized_sources["polint"][nearest_idx] + nearest_fracpol[i] = all_polarized_sources["fracpol"][nearest_idx] + nearest_pol_ra[i] = all_polarized_sources[rm_ra][nearest_idx] + nearest_pol_dec[i] = all_polarized_sources[rm_dec][nearest_idx] # Add columns in delay calibrator table new_columns = { @@ -1090,7 +1167,7 @@ def rmgrid_match( overwrite=True, ) - return polarized_sources + return all_polarized_sources def gaia_quasar_match(file): @@ -1359,6 +1436,7 @@ def generate_catalogues( lotss_catalogue="lotss_catalogue.csv", lbcs_catalogue="lbcs_catalogue.csv", rmgrid_catalogue="../../DR3_RMGrid_v2.fits", + pulsar_catalogue="../../DR3_pulsars_known_v2.fits", pol2delay_radius='0.3deg', lotss_result_file="image_catalogue.csv", delay_cals_file="delay_calibrators.csv", @@ -1604,6 +1682,7 @@ def generate_catalogues( polarized_sources = rmgrid_match( delay_cals_file, rmgrid_catalogue, + pulsar_catalogue, RATar, DECTar, im_radius, @@ -1710,6 +1789,13 @@ def parse_args(): help="input file for rmgrid_catalog [will be downloaded if does not exist]", default="../DR3_RMGrid_v2.fits", # ToDo: Set a default file path/url ) + parser.add_argument( + "--pulsar_catalogue", + dest="pulsar_catalogue", + type=str, + help="input file for pulsar_catalogue [will be downloaded if does not exist]", + default="../DR3_pulsars_known_v2.fits", # ToDo: Set a default file path/url + ) parser.add_argument( "--pol2delay_radius", dest="pol2delay_radius", @@ -1845,6 +1931,7 @@ def main(): lotss_catalogue=args.lotss_catalogue, lbcs_catalogue=args.lbcs_catalogue, rmgrid_catalogue=args.rmgrid_catalogue, + pulsar_catalogue=args.pulsar_catalogue, pol2delay_radius=args.pol2delay_radius, lotss_result_file=args.lotss_result_file, delay_cals_file=args.delay_cals_file, @@ -1906,6 +1993,13 @@ def main(): help="input file for rmgrid_catalog [will be downloaded if does not exist]", default="../DR3_RMGrid_v2.fits", # ToDo: Set a default file path/url ) + parser.add_argument( + "--pulsar_catalogue", + dest="pulsar_catalogue", + type=str, + help="input file for pulsar_catalogue [will be downloaded if does not exist]", + default="../DR3_pulsars_known_v2.fits", # ToDo: Set a default file path/url + ) parser.add_argument( "--pol2delay_radius", dest="pol2delay_radius", @@ -2037,6 +2131,7 @@ def main(): lbcs_catalogue=args.lbcs_catalogue, lotss_result_file=args.lotss_result_file, rmgrid_catalogue=args.rmgrid_catalogue, + pulsar_catalogue=args.pulsar_catalogue, pol2delay_radius=args.pol2delay_radius, delay_cals_file=args.delay_cals_file, match_tolerance=args.match_tolerance, From 154077ff5f89285c9e78e9369d44945f0c201a95 Mon Sep 17 00:00:00 2001 From: HrishikeshShetgaonkar Date: Fri, 28 Aug 2026 09:30:36 +0200 Subject: [PATCH 3/3] Handle different colnames in DR2 DR3 catalogs --- src/plot_field/plot_field.py | 280 +++++++++++++++++++++++------------ 1 file changed, 184 insertions(+), 96 deletions(-) diff --git a/src/plot_field/plot_field.py b/src/plot_field/plot_field.py index 129f2ea..afd6f05 100644 --- a/src/plot_field/plot_field.py +++ b/src/plot_field/plot_field.py @@ -741,38 +741,43 @@ def make_plot( ax.add_patch(c) # plot polarized sources - polint = np.asarray(polarized_sources["polint"], dtype=float) - marker_size = 20 + 100 * polint / np.nanmax(polint) #20 is just a factor to get bigger markers - is_pulsar = polarized_sources["pol_type"] == "Pulsar" - is_rmgrid = ~is_pulsar - - # Normal polarized sources - ax.scatter( - polarized_sources["RA_pol"][is_rmgrid], - polarized_sources["Dec_pol"][is_rmgrid], - transform=ax.get_transform("fk5"), - marker="D", - s=marker_size[is_rmgrid], - facecolors="none", - edgecolors="magenta", - linewidths=1.5, - label="Polarized Source", - zorder=5, - ) + ra = "RA_pol" + dec = "Dec_pol" + if len(polarized_sources) > 0: + polint = np.asarray(polarized_sources["polint"], dtype=float) + marker_size = 20 + 100 * polint / np.nanmax(polint) #20 is just a factor to get bigger markers + is_pulsar = polarized_sources["pol_type"] == "Pulsar" + is_rmgrid = ~is_pulsar + + # LoTSS polarized sources + if len(polarized_sources[ra][is_rmgrid]) > 0: + ax.scatter( + polarized_sources[ra][is_rmgrid], + polarized_sources[dec][is_rmgrid], + transform=ax.get_transform("fk5"), + marker="D", + s=marker_size[is_rmgrid], + facecolors="none", + edgecolors="magenta", + linewidths=1.5, + label="Polarized Source", + zorder=5, + ) - # Pulsars - ax.scatter( - polarized_sources["RA_pol"][is_pulsar], - polarized_sources["Dec_pol"][is_pulsar], - transform=ax.get_transform("fk5"), - marker="D", - s=marker_size[is_pulsar], - facecolors="none", - edgecolors="red", - linewidths=1.5, - label="Pulsars", - zorder=5, - ) + # Pulsars + if len(polarized_sources[ra][is_pulsar]) > 0: + ax.scatter( + polarized_sources[ra][is_pulsar], + polarized_sources[dec][is_pulsar], + transform=ax.get_transform("fk5"), + marker="D", + s=marker_size[is_pulsar], + facecolors="none", + edgecolors="red", + linewidths=1.5, + label="Pulsars", + zorder=5, + ) ### Calculate and plot smallest distance to target dist_ids, dist = smallest_distance(RA, DEC, lbcs) @@ -999,6 +1004,101 @@ def _parse_angular_radius(radius): "Use values such as '1.2deg' or '30arcmin'." ) from exc +def get_rmgrid_catalogue(rmgrid_catalogue=None, outdir="./"): + """ + Return RM-grid catalogue. + If rmgrid_catalogue is supplied, use local file. + Otherwise download the default catalogue and cache it in outdir. + """ + RMGRID_URL = "https://lofar-mksp.org/wp-content/uploads/2022/08/LoTSS_DR2_RMGrid_v1_RMTable.fits_.gz" + # local catalogue supplied + if rmgrid_catalogue is not None: + if os.path.isfile(rmgrid_catalogue): + print( + f"Using user-supplied RM-grid catalogue: " + f"{rmgrid_catalogue}" + ) + return Table.read(rmgrid_catalogue) + else: + raise FileNotFoundError(f"RM-grid catalogue not found: {rmgrid_catalogue}") + + # No catalogue supplied -> check if exist. Else download + rmgrid_catalogue = os.path.join( + outdir, + "LoTSS_DR2_RMGrid_v1_RMTable.fits.gz" + ) + + # Already downloaded? + if os.path.isfile(rmgrid_catalogue): + print( + f"Using cached RM-grid catalogue: " + f"{rmgrid_catalogue}" + ) + + return Table.read(rmgrid_catalogue) + + print("No RM-grid catalogue supplied.") + print(f"Downloading default RM-grid catalogue from {RMGRID_URL}") + + response = requests.get( + RMGRID_URL, + stream=True, + timeout=120, + ) + + response.raise_for_status() + + with open(rmgrid_catalogue, "wb") as f: + for chunk in response.iter_content( + chunk_size=1024 * 1024 + ): + if chunk: + f.write(chunk) + + print( + f"RM-grid catalogue downloaded to: " + f"{rmgrid_catalogue}" + ) + + #Read catalog into memory and then delete the downloaded catalog + rmgrid_table = Table.read(rmgrid_catalogue) + if os.path.isfile(rmgrid_catalogue): + os.remove(rmgrid_catalogue) + + return rmgrid_table + +def standardize_catalogue(table): + """ + Read a table and standardize column names. + Makes life easier when changing catalog versions + """ + aliases = { + "RA_pol": ["RA_pol", "ra_pol", "ra"], + "Dec_pol": ["DEC_pol", "dec_pol", "dec"], + "RM": ["RM", "rm"], + "RM_err": ["RM_err", "rm_err", "rm_error"], + "I_144MHz": ["I_144MHz", "stokesI"], + "polint": ["polint", "PolarizationIntensity"], + "fracpol": ["fracpol", "FractionalPolarization"] + } + + existing = {name.lower(): name for name in table.colnames} + + for standard_name, possible_names in aliases.items(): + for name in possible_names: + if name.lower() in existing: + old_name = existing[name.lower()] + + if old_name != standard_name: + table.rename_column(old_name, standard_name) + + break + + column_names = list(aliases.keys()) + missing = [name for name in column_names if name not in table.colnames] + if missing: + raise KeyError("Missing required column(s) in polarization catalogue: " + ", ".join(missing)) + return table[column_names], column_names def rmgrid_match( delay_cals_file, @@ -1008,29 +1108,19 @@ def rmgrid_match( field_dec, field_radius="1.2deg", pol_radius="0.3deg" -): + ): """find polarized sources within given RA,DEC - first for sources from DR3_RMGrid_v2.fits + first for sources from LoTSS_DR2_RMGrid_v1_RMTable """ - if not os.path.isfile(rmgrid_catalogue): - raise FileNotFoundError(f"RM-grid catalogue not found: {rmgrid_catalogue}") - - rmgrid = Table.read(rmgrid_catalogue) - - rm_ra = "RA_pol" - rm_dec = "Dec_pol" - - pol_columns = ["RA_pol", "Dec_pol", "RM", "RM_err", "polint", "I_144MHz", "fracpol"] - missing = [name for name in pol_columns if name not in rmgrid.colnames] - if missing: - raise KeyError("Missing required RM-grid column(s): " + ", ".join(missing)) + + rmgrid, pol_columns = standardize_catalogue(get_rmgrid_catalogue(rmgrid_catalogue, outdir="./")) radius = _parse_angular_radius(field_radius) field_centre = SkyCoord( float(field_ra), float(field_dec), unit="deg", frame="icrs" ) rm_coords_all = SkyCoord( - rmgrid[rm_ra], rmgrid[rm_dec], unit="deg", frame="icrs" + rmgrid["RA_pol"], rmgrid["Dec_pol"], unit="deg", frame="icrs" ) in_field = rm_coords_all.separation(field_centre) <= radius polarized_sources = rmgrid[in_field] @@ -1041,51 +1131,41 @@ def rmgrid_match( f"({radius.to_value(u.arcmin):.4g} arcmin) of " f"RA={float(field_ra):.6f}, DEC={float(field_dec):.6f}" ) + #Standardize names: + rmgrid_pol = polarized_sources[pol_columns].copy() + rmgrid_pol["pol_type"] = ["RMGrid"] * len(rmgrid_pol) """ now for pulsar sources from DR3_pulsars_known_v2.fits """ - if not os.path.isfile(pulsar_catalogue): - raise FileNotFoundError(f"pulsar catalogue not found: {pulsar_catalogue}") - - pulsars = Table.read(pulsar_catalogue) + pulsar_pol = Table(names=pol_columns) + if pulsar_catalogue is not None and os.path.isfile(pulsar_catalogue): + pulsars, pol_columns = standardize_catalogue(Table.read(pulsar_catalogue)) - pulsar_ra = "ra" - pulsar_dec = "dec" + radius = _parse_angular_radius(field_radius) + field_centre = SkyCoord( + float(field_ra), float(field_dec), unit="deg", frame="icrs" + ) + pulsar_coords_all = SkyCoord( + pulsars["RA_pol"], pulsars["Dec_pol"], unit="deg", frame="icrs" + ) + in_field = pulsar_coords_all.separation(field_centre) <= radius + in_field_pulsars = pulsars[in_field] - pol_columns = ["ra", "dec", "RM", "RM_err", "polint", "I_144MHz", "fracpol"] - missing = [name for name in pol_columns if name not in pulsars.colnames] - if missing: - raise KeyError("Missing required pulsar column(s): " + ", ".join(missing)) + print( + f"Found {len(in_field_pulsars)} pulsars " + f"{radius.to_value(u.deg):.4g} deg " + f"({radius.to_value(u.arcmin):.4g} arcmin) of " + f"RA={float(field_ra):.6f}, DEC={float(field_dec):.6f}" + ) - radius = _parse_angular_radius(field_radius) - field_centre = SkyCoord( - float(field_ra), float(field_dec), unit="deg", frame="icrs" - ) - pulsar_coords_all = SkyCoord( - pulsars[pulsar_ra], pulsars[pulsar_dec], unit="deg", frame="icrs" - ) - in_field = pulsar_coords_all.separation(field_centre) <= radius - in_field_pulsars = pulsars[in_field] + #Standardize names: + pulsar_pol["pol_type"] = ["Pulsar"] * len(pulsar_pol) - print( - f"Found {len(in_field_pulsars)} pulsars " - f"{radius.to_value(u.deg):.4g} deg " - f"({radius.to_value(u.arcmin):.4g} arcmin) of " - f"RA={float(field_ra):.6f}, DEC={float(field_dec):.6f}" - ) + else: + print(f"pulsar catalogue not found: {pulsar_catalogue}. Continuing without") # combine both catalogs - rmgrid_pol = polarized_sources[["RA_pol", "Dec_pol", "RM", "RM_err", "polint", "I_144MHz", "fracpol"]].copy() - rmgrid_pol["pol_type"] = ["RMGrid"] * len(rmgrid_pol) - # Pulsars - pulsar_pol = in_field_pulsars[ - ["ra", "dec", "RM", "RM_err", "polint", "I_144MHz", "fracpol"] - ].copy() - - pulsar_pol.rename_column("ra", "RA_pol") - pulsar_pol.rename_column("dec", "Dec_pol") - pulsar_pol["pol_type"] = ["Pulsar"] * len(pulsar_pol) # Combine them if len(rmgrid_pol) > 0 and len(pulsar_pol) > 0: @@ -1114,8 +1194,8 @@ def rmgrid_match( # find polarized sources near the delay cals if len(all_polarized_sources) > 0: pol_coords = SkyCoord( - all_polarized_sources[rm_ra], - all_polarized_sources[rm_dec], + all_polarized_sources["RA_pol"], + all_polarized_sources["Dec_pol"], unit="deg", frame="icrs", ) @@ -1142,8 +1222,8 @@ def rmgrid_match( nearest_pol_dist[i] = separations[nearest_idx].deg nearest_polint[i] = all_polarized_sources["polint"][nearest_idx] nearest_fracpol[i] = all_polarized_sources["fracpol"][nearest_idx] - nearest_pol_ra[i] = all_polarized_sources[rm_ra][nearest_idx] - nearest_pol_dec[i] = all_polarized_sources[rm_dec][nearest_idx] + nearest_pol_ra[i] = all_polarized_sources["RA_pol"][nearest_idx] + nearest_pol_dec[i] = all_polarized_sources["Dec_pol"][nearest_idx] # Add columns in delay calibrator table new_columns = { @@ -1435,8 +1515,8 @@ def generate_catalogues( bright_limit_Jy=5.0, lotss_catalogue="lotss_catalogue.csv", lbcs_catalogue="lbcs_catalogue.csv", - rmgrid_catalogue="../../DR3_RMGrid_v2.fits", - pulsar_catalogue="../../DR3_pulsars_known_v2.fits", + rmgrid_catalogue=None, + pulsar_catalogue=None, pol2delay_radius='0.3deg', lotss_result_file="image_catalogue.csv", delay_cals_file="delay_calibrators.csv", @@ -1708,6 +1788,8 @@ def generate_catalogues( fit_spectrum(delay_cals_file, outdir) if vlass: + from plot_field.vlass_search import search_vlass + ## Get cutouts of all LBCS sources print("Getting cutouts of LBCS sources") for i, source in enumerate(result): @@ -1726,6 +1808,10 @@ def generate_catalogues( convert_vlass_fits(outfile) except Exception: print(f"VLASS Download failed for {source['Observation']}") + search_vlass(c, crop=True, crop_scale=256) + os.system("mv vlass_post**.fits %s" % outfile) + convert_vlass_fits(outfile) + if html: app = make_html( RATar, @@ -1741,6 +1827,8 @@ def generate_catalogues( ) app.run_server(debug=True, use_reloader=False) + # return + def parse_args(): parser = argparse.ArgumentParser( @@ -1786,22 +1874,22 @@ def parse_args(): "--rmgrid_catalogue", dest="rmgrid_catalogue", type=str, - help="input file for rmgrid_catalog [will be downloaded if does not exist]", - default="../DR3_RMGrid_v2.fits", # ToDo: Set a default file path/url + default=None, + help="input file for rmgrid_catalog [will be downloaded if does not exist]" ) parser.add_argument( "--pulsar_catalogue", dest="pulsar_catalogue", type=str, - help="input file for pulsar_catalogue [will be downloaded if does not exist]", - default="../DR3_pulsars_known_v2.fits", # ToDo: Set a default file path/url + default=None, + help="input file for pulsar_catalogue [will be downloaded if does not exist]" ) parser.add_argument( "--pol2delay_radius", dest="pol2delay_radius", type=str, help=( - "Radius around delay calibrator to choose polarized sources. " + "Radius around delay calibrators to choose polarized sources. " "Examples: 0.3deg, 30arcmin [default: 0.3deg]" ), default="0.3deg", @@ -1990,15 +2078,15 @@ def main(): "--rmgrid_catalogue", dest="rmgrid_catalogue", type=str, - help="input file for rmgrid_catalog [will be downloaded if does not exist]", - default="../DR3_RMGrid_v2.fits", # ToDo: Set a default file path/url + default=None, + help="input file for rmgrid_catalog [will be downloaded if does not exist]" ) parser.add_argument( "--pulsar_catalogue", dest="pulsar_catalogue", type=str, - help="input file for pulsar_catalogue [will be downloaded if does not exist]", - default="../DR3_pulsars_known_v2.fits", # ToDo: Set a default file path/url + default=None, + help="input file for pulsar_catalogue [will be downloaded if does not exist]" ) parser.add_argument( "--pol2delay_radius",