Versions Compared

Key

  • This line was added.
  • This line was removed.
  • Formatting was changed.

...

  1. 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).
  2. 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.
  3. 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.
    1. 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 actually stores 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

Here's a function from J-drive in ".../CIRA/Operations/Maps&Locations/MWA/MWA_Ph2_Extended-LB/Long Baselines Coordinate conversions":

...

languagepy

...

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":

Code Block
languagepy
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, 64) / 72024

    lof1fact4 = _a1(61 /+ (n090 * k0)t0 + 298 * Q0 lof2+ =45 (1 + 2* t0 * t0 +- Q0)252 * math.pow(dd0, 3) / 6.0e1sq - 3 * Q0 * Q0) * math.pow(dd0, 6) / 720

    lof3lof1 = _a1 / (5 - 2 * Q0n0 * k0)
    lof2 = (1 + 282 * t0 -+ 3Q0) * math.pow(Q0, 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 surveyed coords to lat, lon and height

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 attached. 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

...

 * zone - 183.0) or 3.0) - _a3

    return (latitude, longitude)


Converting from lat, lon and height to X,Y,Z

...

  1. 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

  2. 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
  3. 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 actually stores 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. 

...