uvw errors for non-terrestrial telescopes?
Closed this issue · 6 comments
In our telecon the other day, we looked at the calculation of uvw's (which eventually go into the visibility calculation) for telescopes on the moon. When initializing from a parameter file, the UVData method set_uvws_from_antenna_positions gets called, which I thought we decided calls code from utils.pyx (i.e. cython code) where gps parameters about the Earth curvature and such are hardcoded. Taking another look, neither @aelanman nor I couldn't find where that method called the cython code. @mkolopanis Can you remember where the issue was? Am I misremembering which method called the cython with hardcoded earth parameters?
functions which depend on the hardcoded GPS constants are the _lla_to/from_xyz conversions
https://github.com/RadioAstronomySoftwareGroup/pyuvdata/blob/119315929ed698ac735bacd868612c082be61a23/pyuvdata/utils.pyx#L213
so anything in pyuvdata that needs LatLonAlt_from_XYZ or XYZ_from_LatLonAlt depends on them. which includes things like get_enu_antpos and the ECEF -> ENU conversion functions
Ok, so maybe I was wrong that it was the set_uvws_from_antenna_positions method in particular. But it looks like parse_telescope_params on line 1220 of simsetup.py runs:
tele_params['telescope_location'] = uvutils.XYZ_from_LatLonAlt(*telescope_location)
So that is likely an issue. Not sure if it's the only one...
Oh, that's subtle... pyuvsim.simsetup sets the XYZ assuming the Lat/Lon/Alt are on Earth, and then SkyModel uses the XYZ position as selenocentric coordinates.
A quick way to check if that's what's happening is to look at the telescope_location in the output file and see if the position vector is roughly an Earth radius vs. a Moon radius.
I can't say for certain, but my guess is the output location is an Earth radius like location. I say that because I believe pyuvdata has a check that the radius is "sensible" and we don't fail that check. But we can look to be sure.
After an initial conversation with @aelanman this morning, we decided there are two issues, one easier to fix than the other:
(1) line 1220 in simsetup sets the telescope location (that will be attached to the UVData object) as ECEF even when presented with a lunar lat/lon/height. To fix, we should check the world and if it's the Moon, call the MoonLocation.from_selenocentric method to set the XYZ coordinates on the UVData object. (Whether we want to call the method from the object or turn the last couple lines of code in that method into a standalone function in LunarSky is debatable.) This fix will also likely trigger a check in pyuvdata about the distance from the center of the earth to fail.
(2) line 1258 in simsetup calculates the antenna positions in ECEF. It then subtracts off the telescope location, which may offset this error. However, once issue (1) is fixed, this will definitely be an error. To fix, we need to calculate the telescope locations used the LunarTopo transform lunartopo_to_mcmf. As written, LunarTopo objects are initialized in an alt/az like coordinate system, but in principle can be represented in an ENU coordinate system. @aelanman may need to make some changes to LunarSky to enable this transformation to proceed as we need it to.
I've started working to address this (see branch fixing_lunar_pos). Making the replacement to correct (1) did cause errors in the pyuvdata acceptability checks. I've disabled the check for now. I made a change to fix (2) as well, but that's causing a failure in pyuvdata because it's using ENU_from_ECEF to reverse the transformation I'm doing here (and so not liking that the vectors aren't Earth-radius in scale). It may be possible to fix that using utils.py:calc_uvws directly.
Overall, it seems like these issues can be resolved quickly with some patches to pyuvsim, but long-term it would make more sense to properly incorporate support for MCMF and LunarTopo frames in pyuvdata directly using the telescope_frame parameter (as I think was the idea in the pyuvdata slack).