WCS Coordinates¶

WCS (World Coordinate System) information is preserved during FITS-to-Zarr conversion, enabling celestial coordinate (RA/Dec) overlays on sky images. The radport accessor provides methods to check WCS availability, convert between pixel and sky coordinates, and create publication-quality WCS-projected plots.

Checking WCS Availability¶

import ovro_lwa_portal as ovro

ds = ovro.open_dataset("path/to/data.zarr")

if ds.radport.has_wcs:
    print("WCS coordinates available")

The WCS header is searched in three locations (in order):

Variable attributes (fits_wcs_header on the SKY variable)
Dataset attributes (fits_wcs_header on the dataset)
A wcs_header_str data variable in the dataset

!!! note WCS functionality requires astropy to be installed.

Pixel-to-Sky Conversion¶

Convert pixel indices along the l and m dimensions to celestial coordinates at a specific observation time. Pass exactly one of time_idx (index into the dataset time coordinate) or time_mjd (MJD; nearest time is used):

ra, dec = ds.radport.pixel_to_coords(l_idx=100, m_idx=100, time_idx=0)
# or: ra, dec = ds.radport.pixel_to_coords(100, 100, time_mjd=60000.0)
print(f"RA: {ra:.4f} deg, Dec: {dec:.4f} deg")

RA is returned in the range [0, 360) degrees
Dec is returned in degrees
Requires a time dimension on the dataset
Raises ValueError if indices are out of bounds, if both or neither of time_idx and time_mjd are given, or if the direction is below the horizon at that time

Sky-to-Pixel Conversion¶

Convert RA/Dec coordinates to pixel indices:

l_idx, m_idx = ds.radport.coords_to_pixel(ra=180.0, dec=45.0, time_idx=0)
print(f"Pixel: l={l_idx}, m={m_idx}")

Returns rounded integer indices
Raises ValueError if the coordinates fall outside the image bounds
Uses the FK5 celestial frame

Plotting with WCS Projection¶

The plot_wcs method creates plots with RA/Dec coordinate axes and a grid overlay:

# Basic WCS plot
fig = ds.radport.plot_wcs(freq_mhz=50.0)

Customizing the Plot¶

# Publication-quality dark-background plot
fig = ds.radport.plot_wcs(
    freq_mhz=50.0,
    mask_radius=1800,
    cmap="inferno",
    grid_color="white",
    grid_alpha=0.6,
    grid_linestyle=":",
    label_color="white",
    facecolor="black",
    figsize=(10, 10),
)

Key Parameters¶

Parameter	Default	Description
`var`	`"SKY"`	Data variable to plot (`"SKY"` or `"BEAM"`)
`time_idx`	`0`	Time index
`freq_idx`	`0`	Frequency index (ignored if `freq_mhz` is set)
`freq_mhz`	`None`	Select frequency by value in MHz
`cmap`	`"inferno"`	Matplotlib colormap
`robust`	`True`	Use 2nd/98th percentile for color scaling
`mask_radius`	`None`	Circular mask radius in pixels
`grid_color`	`"white"`	Color of the RA/Dec grid lines
`grid_alpha`	`0.6`	Transparency of grid lines
`facecolor`	`"black"`	Plot background color
`add_colorbar`	`True`	Whether to include a colorbar

The plot automatically handles RA axis inversion (RA increases to the left on the sky) and uses robust percentile-based scaling by default.

Example Workflows¶

Locate a Source by RA/Dec and Extract a Light Curve¶

# Find the Crab Nebula (RA=83.633, Dec=22.014)
l_idx, m_idx = ds.radport.coords_to_pixel(ra=83.633, dec=22.014, time_idx=0)

# Get the l, m coordinate values at those indices
l_val = float(ds.coords["l"].values[l_idx])
m_val = float(ds.coords["m"].values[m_idx])

# Extract the light curve as a DataArray, or plot with the accessor (same l, m)
lc = ds.radport.light_curve(l=l_val, m=m_val)
fig = ds.radport.plot_light_curve(l=l_val, m=m_val)

Create a Cutout Around a Known Position¶

# Convert RA/Dec to pixel coordinates
l_idx, m_idx = ds.radport.coords_to_pixel(ra=83.633, dec=22.014, time_idx=0)
l_val = float(ds.coords["l"].values[l_idx])
m_val = float(ds.coords["m"].values[m_idx])

# Extract a cutout centered on the source
cutout = ds.radport.cutout(l_center=l_val, m_center=m_val, dl=0.1, dm=0.1)
fig = ds.radport.plot_cutout(cutout)

WCS Plot with Detected Sources¶

import matplotlib.pyplot as plt

# Create WCS plot
fig = ds.radport.plot_wcs(freq_mhz=50.0, mask_radius=1800)
ax = fig.axes[0]

# Overlay detected sources
peaks = ds.radport.find_peaks(time_idx=0, freq_idx=0, threshold=5.0)
for peak in peaks:
    ra, dec = ds.radport.pixel_to_coords(
        int(peak["l_idx"]), int(peak["m_idx"]), time_idx=0
    )
    ax.plot(ra, dec, "r+", markersize=10, transform=ax.get_transform("fk5"))

plt.show()

How WCS Is Preserved¶

The WCS pipeline works as follows:

Original FITS files contain standard WCS headers (CRVAL, CRPIX, CDELT, etc.)
Header fixing (fix_fits_headers) corrects any non-standard headers
Zarr conversion stores the complete WCS header as a string attribute (fits_wcs_header) on the data variable or dataset
At load time, the radport accessor parses the stored header string back into an astropy WCS object on demand

This means WCS information survives the FITS-to-Zarr conversion pipeline and is available whenever the original FITS files contained valid WCS headers.