Fix FITS headers for improved astrometry

[dentalfloss1]

The wcs objects that we have been using, in particular in the FITS header, is almost correct, but generates positions that are off by typically about a pixel. Logan Cordonnier has a working version generated from fitting positions. It would be nice to have a rigorous theoretical derivation of the wcs object, but in the interest of time it may be better just to incorporate what he has written.

@jaycedowell let me know if you have an opinion on this matter.

[jaycedowell]

For now we can simply add Logan corrections if we have an appropriate mechanism to only apply them to the LWA-SV data. Ultimately, though, we should understand how to take a shifted phase center and convert that into the correct WCS parameters.

[jaycedowell]

@lcordonnier what all needs to change in the WCS for Sevilleta data?

[lcordonnier]

I've optimized 5 parameters for Sevilleta:

• Lonpole
• Center HA and Dec (since HA stays fixed)
• xi and eta (to control projection center, given by keywords PV2_1a, PV2_2a)

The best values I got for these, respectively:

• 179.21441725378727
• 357.38856977271047, 33.507121493107995
• 0.011830877131833636, -0.008194262566706176

Granted, my sky image was flipped for this (NESW vs standard NWSE, going clockwise) though I don't think that should affect the values.

[jaycedowell]

orville_imager.py uses:

self.phase_center_ha = 1.0*ephem.hours("-0:07:59.82")
self.phase_center_dec = 1.0*ephem.degrees("33:21:27.5")

and the station is taken to be at:

GEO_N   +34.348358 # N Lat [deg]
GEO_E  -106.885783 # E Lat [deg]
GEO_EL 1477.8   #  [m] above MSL - FIXME

[jaycedowell]

I found my old code called optimalPointing.py that's used to find the optimal phase center of a station. That HA/δ combo corresponds to moving the phase center by 1.933° at in the direction of an azimuth of 120.65° at LWA-SV.

[jaycedowell]

@lcordonnier I was looking at Calabretta & Greisen (2002), Section 5.1.5 again and I noticed that there looks to be a dependency of ξ and η on the center of the image (θ and φ). I think that θ and φ should be related to the azimuth/elevation of the phase center which, in turn, is related to the phase center HA and δ. When you fit did you treat all of the parameters are independent?

[lcordonnier]

I let it optimize all five parameters independently, yes. So it didn't explicitly account for a dependency between ξ & η and center HA & δ, though if the relationship between them is known I can go back and refit taking that into account. In actuality, the θ and φ parameters were directly used in the optimization, which were then converted to ξ and η for the WCS using Eqns. 63/64 in the linked paper.

[jaycedowell]

Maybe try a fit where where the only free parameters are θ, φ, and LONGPOLE. We can then find HA_c, δ_c, ξ, and η from those.

The only thing that I don't really know where it fits in is LONGPOLE. It's got to be related to something. How sensitive are the fits to that parameter?

[lcordonnier]

The lonpole value has ranged from 179.15 - 179.48 across the ~10 different fits, which span between 1-4 days of data and 3-7 free parameters.

I can run the fit again with only those 3 free parameters, however I'll need to use some value for HA_c and δ_c in order to generate the WCS; should I assume that θ = alt_c, φ = az_c -> convert using astropy -> HA_c, δ_c or is the relationship more complex than that?

[jaycedowell]

I can run the fit again with only those 3 free parameters, however I'll need to use some value for HAc and δc in order to generate the WCS; should I assume that θ = alt_c, φ = az_c -> convert using astropy -> HA_c, δ_c or is the relationship more complex than that?

That looks reasonable. We may also want to make sure EQUINOX/EPOCH is set correctly for the WCS. They may need to be epoch-of-date rather than J2000.

[jaycedowell]

Random thought: We could take a TBF capture with the correct frequency range, run it through correlateTBF.py to generate a measurement set, and then image that with wsclean/the wsclean (u,v) shifter to see what that WCS looks like.

[lcordonnier]

Ok, so I've tried tricking WCS into using Az/Alt coordinates in place of RA/Dec. This has a couple of benefits:

1. EQUINOX/EPOCH shouldn't be a concern since the Az/Alt coordinates natively take into account the time of observation (I assume that WCS doesn't try to precess things if EQUINOX/EPOCH aren't explicitly given, but I'm not 100% sure about this. Since the coordinates are still technically RA/Dec in its eyes, it may be doing something to the coordinates under the hood that I'm not aware of.)
2. I can now treat θ = alt_c and φ = az_c directly without having to convert them to HA/Dec with astropy, again avoiding J2000 vs epoch-of-date concerns.

I've optimized the positions with just LONPOLE, az_c, and alt_c and found the best-fit parameters to be, respectively:

209.6121637218684, 118.76647246379336, 87.5850923296785

Where the corresponding header values would be:

180.0, 90.97437629624612, 88.06662829704536

Perhaps it's just a coincidence, but converting the phase center HA/Dec value from orville_imager.py (assuming it's measured in epoch-of-date) into Az/Alt gives a coordinate of:

120.26495457, 88.0666283

which shows much better agreement with the best-fit Az coord. Again, may just be a coincidence. The offset in the best-fit LONPOLE from the expected 180 (~29.6) is similar to the offset in best-fit az_c from header az_c (~27.8). I'd speculate that these offsets being the same (or close to the same) has to do with using Az/Alt coords in WCS (since now the 'celestial pole' that LONPOLE is measured w.r.t. is actually the zenith). This may also just be a coincidence though.

The quality of the optimization is... okay. The pixel offset rms is 2x higher than it is for optimizing all 5 parameters independently (0.45 px vs 0.24 px), but half that (0.45 px vs 0.81 px) of just using the values directly from the header (i.e. before I started optimizing the astrometry). Here's an image showing the comparison between my best astrometry so far (5 parameter) and the best-fit using the Az/Alt method:

[jaycedowell]

That's an interesting result that the HA/dec coordinates do not convert to the same values for az/alt that are in orville_imager.py. How did you do that conversion?

[lcordonnier]

Here's the code snippet I used to do the conversion:

from astropy.coordinates import Angle
import astropy.coordinates as ac
from astropy.time import Time
import astropy.units as u

HA = Angle('-0h07m59.82s')
Dec = Angle('33d21m27.5s')

SV = ac.EarthLocation(lat=34.348358*u.deg, lon=-106.885783*u.deg, height=1477.8*u.m)
t_obs = Time(59670.583333333336, format='mjd') #time of observation
altaz_obj = ac.AltAz(location=SV, obstime=t_obs)

# assume the HA/Dec are measured in the epoch-of-date
hadec = ac.SkyCoord(ha=HA, dec=Dec, frame='hadec', obstime=t_obs, location=SV)
altaz = hadec.transform_to(altaz_obj)
print(altaz)

Which outputs:

<SkyCoord (AltAz: obstime=59670.583333333336, location=(-1531554.7717322, -5045440.98395596, 3579249.98860634) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt) in deg
    (120.26495457, 88.0666283)>