MWA coordinates and conversion
Radio telescopes require a well defined coordinate systems for core calculations and to make accurate maps of the sky. The three core systems that have to be well defined are the celestial system (and epoch), the time system, and the earth-based system describing the location of the array antennas. Typically an array centre is defined, which defines the origin of a local cartesian coordinate system, and defines the latitude and longitude of the array.
MWA tiles were surveyed by professional surveyors, who present the results in UTM coordinates "easting", "northing" and "height". This is a standard system for surveyors, but the directions for "easting" and "northing" don't correspond to local east and north in general. We need to convert the surveyed coords to a standard cartesian coordinate system used by radio telescopes.
It is important to note that the 128T original tiles were surveyed at the south-west corner of the tile, so the true centre of the tile is +2.5 m to the north and east of this point. This procedure was also used for the Long Baseline (solar powered) tiles, however the "Hexes" may not have done this, so there are corrections that need to me made to the raw surveyed coordinates for all to be consistent.
The basic recipe for conversion is as follows:
- Convert the surveyed coords to latitude, longitude and height. This system needs a reference ellipsoid such as the WGS84 system, and the height is assumed to be relative to the reference ellipsoid (it is hard to find definitive information on what the height actually means).
- Convert the lat, lon and height to the X,Y,Z absolute "Earth centred, Earth fixed" coordinate system where the origin is the centre of the earth. This conversion needs to assume the same ellipsoid as the conversion from UTM coords. The X axis points out of the equator at zero degs longitude, the Z axis points through the north pole along the rotation axis of the earth, and the Y axis completes a right-handed coordinate system. (Y points out of the equator in the Indian ocean at longitude 90 degs.) It is possible to specify the location of the telescope antennas in this X,Y,Z system in many software packages, but the numbers are very large and non-intuitive, so a local coordinate system is often preferred.
- Convert the X,Y,Z coordinates of the antennas to a local cartesian coordinate system where the origin of the system is the defined centre of the array. The axes of the local system are east (defined at the array centre), north (defined at the array centre) and height or "up". It is also worth noting that some software confusingly defines another "local X,Y,Z" system where the global X,Y,Z system has been rotated around the Z axis such that X points out of the equator at the longitude of the array centre. I'm not sure if there is an official name to this system, and it doesn't help much because the coordinates of the antennas are still large numbers and non-intuitive because the local horizon is angled to all coordinate axes.
- To convert to local E,N,U units, a rotation on the sphere is required. Note that rotations on sphere do not commute, and we need to do the longitude rotation before the latitude rotation when going from X,Y,Z to E,N,U. If going from E,N,U to X,Y,Z, then the opposite order would be used.
Note that the MWA tile location database seems to actually store E, N and the raw height above the reference ellipsoid (the original survey height value), NOT the height above the tile center (calculated as 'U' above).
Converting surveyed UTM coordinates to lat/lon
Converting from UTM coords to lat, lon and height uses well established formuae like Redfearn's formula. There are online tools that allow code to be checked against trusted sources. Example code is below. Just reiterating, the height is assumed to be height above the reference ellipsoid.
E.g. for MWA centre in UTM Map Grid Australia (MGA) zone 50:
cent_easting = 467254.490961539 # meters
cent_northing= 7046381.90073077 # meters
cent_alt = 377.8269 # meters
The array centre is: lat -26.70331940, lon 116.67081524 degs.
Here's a function from J-drive in ".../CIRA/Operations/Maps&Locations/MWA/MWA_Ph2_Extended-LB/Long Baselines Coordinate conversions":
def utmToLatLng(zone, easting, northing, northernHemisphere=True):
if not northernHemisphere:
northing = 10000000 - northing
a = 6378137
e = 0.081819191
e1sq = 0.006739497
k0 = 0.9996
arc = northing / k0
mu = arc / (a * (1 - math.pow(e, 2) / 4.0 - 3 * math.pow(e, 4) / 64.0 - 5 * math.pow(e, 6) / 256.0))
ei = (1 - math.pow((1 - e * e), (1 / 2.0))) / (1 + math.pow((1 - e * e), (1 / 2.0)))
ca = 3 * ei / 2 - 27 * math.pow(ei, 3) / 32.0
cb = 21 * math.pow(ei, 2) / 16 - 55 * math.pow(ei, 4) / 32
cc = 151 * math.pow(ei, 3) / 96
cd = 1097 * math.pow(ei, 4) / 512
phi1 = mu + ca * math.sin(2 * mu) + cb * math.sin(4 * mu) + cc * math.sin(6 * mu) + cd * math.sin(8 * mu)
n0 = a / math.pow((1 - math.pow((e * math.sin(phi1)), 2)), (1 / 2.0))
r0 = a * (1 - e * e) / math.pow((1 - math.pow((e * math.sin(phi1)), 2)), (3 / 2.0))
fact1 = n0 * math.tan(phi1) / r0
_a1 = 500000 - easting
dd0 = _a1 / (n0 * k0)
fact2 = dd0 * dd0 / 2
t0 = math.pow(math.tan(phi1), 2)
Q0 = e1sq * math.pow(math.cos(phi1), 2)
fact3 = (5 + 3 * t0 + 10 * Q0 - 4 * Q0 * Q0 - 9 * e1sq) * math.pow(dd0, 4) / 24
fact4 = (61 + 90 * t0 + 298 * Q0 + 45 * t0 * t0 - 252 * e1sq - 3 * Q0 * Q0) * math.pow(dd0, 6) / 720
lof1 = _a1 / (n0 * k0)
lof2 = (1 + 2 * t0 + Q0) * math.pow(dd0, 3) / 6.0
lof3 = (5 - 2 * Q0 + 28 * t0 - 3 * math.pow(Q0, 2) + 8 * e1sq + 24 * math.pow(t0, 2)) * math.pow(dd0, 5) / 120
_a2 = (lof1 - lof2 + lof3) / math.cos(phi1)
_a3 = _a2 * 180 / math.pi
latitude = 180 * (phi1 - fact1 * (fact2 + fact3 + fact4)) / math.pi
if not northernHemisphere:
latitude = -latitude
longitude = ((zone > 0) and (6 * zone - 183.0) or 3.0) - _a3
return (latitude, longitude)
Converting from lat, lon and height to X,Y,Z
This is analogous to finding the cartesian coords of a point on a sphere, except it is on the ellipsoid so there are some additional terms. Example python code based on WGS84 ellipsoid:
def Geodetic2XYZ(lat_rad,lon_rad,height_meters):
EARTH_RAD_WGS84 = 6378137.0 # meters in the WGS84 (effectively same as GRS80)
E_SQUARED = 6.69437999014e-3
s_lat = math.sin(lat_rad)
c_lat = math.cos(lat_rad)
s_lon = math.sin(lon_rad)
c_lon = math.cos(lon_rad)
chi = math.sqrt(1.0 - E_SQUARED*s_lat*s_lat)
X = (EARTH_RAD_WGS84/chi + height_meters)*c_lat*c_lon
Y = (EARTH_RAD_WGS84/chi + height_meters)*c_lat*s_lon
Z = (EARTH_RAD_WGS84*(1.0-E_SQUARED)/chi + height_meters)*s_lat
return (X,Y,Z)
This calculation needs to be done for all antennas, and the centre of the array.
E.g. the centre of the MWA in these coords for Array centre: lat -26.70331940, lon 116.67081524 is
X: -2559454.08, Y: 5095372.14, Z: -2849057.18
Converting from X,Y,Z to E,N,U
We've defined an array centre. During the early days of MWA, a surveying mark was made on a rock outside the MWA donga, and this point was a well-defined and accurately measured reference point for the surveyors. Since this point is close to the hub of the MWA, it is as good a point as any to define as the centre of the array, so it was.
- Calculate the coords of the antennas in X,Y,Z relative to the array centre's X,Y,Z: this is just a vector difference:
ant_geo_X = ant_X - arr_X
ant_geo_Y = ant_Y - arr_Y
ant_geo_Z = ant_Z - arr_Z - Rotate the longitude:
x_local = x_geo*math.cos(lon_rad) + math.sin(lon_rad)*y_geo
y_local = y_geo*math.cos(lon_rad) - math.sin(lon_rad)*x_geo
z_local = z_geo - Rotate the latitude:
e_local = y_local
n_local = z_local*math.cos(lat_rad) - math.sin(lat_rad)*x_local
u_local = x_local*math.cos(lat_rad) + math.sin(lat_rad)*z_local
Note that the MWA tile location database seems to actually store E, N and the raw height above the reference ellipsoid (the original survey height value), NOT the height above the tile center (calculated as 'u_local' above).
Some versions of these codes have been put on the J drive under "EECMS/CIRA/Operations/Maps&Locations/MWA", but these are clearly copies of Randall's original code that he had no idea were there, so they should not be considered authoritative.
History
32T
Tile locations for the prototype 32T array were surveyed sometime in 2008. The most definitive document available for the locations is "32T_Reordered_Tile_Corner_Coordinates.xls" in the Curtin J drive under .../MWA/Phase0-32T/Construction. It isn't known if this is truly the original, or derived from something else. This includes the location of the marker rock outside the MWA donga (these match exactly the numbers for the MWA centre above). The spreadsheet doesn't say who did the surveying or what the coordinate system was, but putting the easting and northing values into a Geoscience Australia Converter gives the exact public coordinates for the WGS84 and GRS80 ellipsoids, so it is reasonable to assume WGS84 was used.
Phase I tiles
The original survey data spreadsheet for the tiles appears to be on the J drive in .../CIRA/Operations/MWA/Phase1-128T/Locations/
The file "128_Surveyed_Locations.xls" appears to be the original from HTD, the surveying company. It does not say so explicitly in this spreadsheet, but the surveyed points were the south-west corners of the tiles, so the true centre of the tile is +2.5m north and east of that point. The surveyed locations of the receivers (or receiver pads) is in a file "ASCON Rx Locations.xls" in the same directory.
The extracted tileID, Easting, northing and RL extracted as text from the spreadsheet are in the file "ant_pos_128T_UTM.txt" in .../CIRA/Operations/Maps&Locations/MWA/MWA_Ph1_128T/Conversions on the J drive. The python2 script "conv_128t_ORIGINAL.py" in this same directory was used to convert the coords in ant_pos_128T_UTM.txt to XYZ and ENU coords.
Dave Emrich notes regarding the original survey data (16 Sep 2022): There is the possibility that one of the tiles on the east side of the “core” was “flipped” after this survey was completed, so that the HTD survey point in this file, actually reflects the SE corner, but that should be easily identified by comparing the E/N values in this spreadsheet with the ‘current’ ones, … an error of 5m true-East should pop out!
Phase II hex tiles
From Dave Emrich:
"So far as I remember the hex tiles only had five coordinates professionally surveyed, a centre, and four cardinals 50m? or so true N,S,E,W from the centre. One of the many spreadsheets should contain the HTD survey information for those ten co-oordinate points (five for each hex). It may also have the flag-pole and other locations in it. (I feel sure we had HTD do more than just the 10 points for the hexes, but I can't remember what). The locations of the individual hex tiles were computed as offsets from those cardinals and set out "manually" by us. Given the entire space was cleared of any vegetation and large rocks, there wasn't any need to move away from the 'desired' locations, so there wasn't a need to survey them after they were placed. (There was some initial mucking around with Hex S1 due to a slight rise in ground level but that was subsequently evened out and the tile set in its proper "desired" location)."
There's one file ("HexesRev4_TST_Output_Edited.txt" in ".../Maps&Locations/MWA/MWA_Ph2_Compact-Hexes/Hex Coordinate conversions/rev4" that has UTM values that match (to within a few cm) the values in the database, after adding 2.5m to E and N values (to convert SE corner to tile centers).
Phase II long baseline tiles
The raw survey data appears to be in the spreadsheet "HTD LONG BASELINE SURVEY.XLSX" in .../CIRA/Operations/Maps&Locations/MWA/Surveys/, this is from the surveying company HTD. The spreadsheet says "SW Corner POINTS DATA", so this is again the location of the SW corner of the tile. This spreadsheet also notes the coord convention: MGA94 Z50.
This file appears to be duplicated as "MWA_PH2-LB_Surveyed_Coordinates-HTD.xlsx" in .../CIRA/Operations/Maps&Locations/MWA/MWA_Ph2_Extended-LB/ and the values extracted out of it into text as the file "LB_Survey_Output.txt".