CODEWITHSUNDEEP
pythonmap-projections
Drakekin
Drakekin

Reputation: 1348

Converting a point on a TOAST map to the lat/lon on a sphere

I have an image of a map of a planet in the TOAST projection of size n, where n is a power of two. I would like to take a pixel on the TOAST map and calculate the latitude/longitude of the corresponding point on a sphere.

I have the following function, but it does not appear to work.

def _x_y_to_lat_lon((x, y), size):
    quad = size / 2

    if x >= quad:
        x -= quad
        if y >= quad:
            quadrant = 1
            y -= quad
        else:
            quadrant = 2
            y = (quad - 1 - y)
    else:
        x = (quad - 1 - x)
        if y >= quad:
            quadrant = 0
            y -= quad
        else:
            quadrant = 3
            y = (quad - 1 - y)

    x += 1
    y += 1

    ax = abs(x)
    ay = abs(y)
    if ax + ay > size:
        hemisphere = -1
    else:
        hemisphere = 1

    latitude = (ax + ay) / size * pi * hemisphere
    longitude = (ax - ay) / (ax + ay) * pi / 2 * hemisphere + (pi / 4) + (quadrant * pi / 2)
    return latitude, longitude

Upvotes: 0

Views: 253

Answers (2)

han.k
han.k

Reputation: 21

The second solution is also not fully correct. It gives a strange non linear longitude. Longitudes around 45, 135, 225 and 315 degrees are too much expanded and at 0,90,180 and 270 degrees longitude are too much shrinked. Reason is the assumption that the longitude in the TOAST projection is the angle equals atan2(u/v) This is not correct. The longitude a the distance of the line with constant latitude, not an angle. See the sketch.

If the origin (0,0) of the TOAST projection is in the middle and X:=V and Y:=V and are both positive, the RA, DEC for the star in the sketch is then:

latitude or sky declination:=(y+x-1)*pi/2;

longitude or sky right_ascension:=0.5*pi*(1-x/(y+x));

Apply this formula for all 8 triangles in the TOAST projection with each time the correct correction as follows:

(1) First make X, Y positive and in range [0..1] from the pole.

(2) Apply above formula

(3) Apply Lon/DEC and lat/RA correction/shift each time for each triangle.

Noted that the Microsoft TOAST program doesn't make straight declination/latitude lines near the poles. My routine does.

TOAST projection of the sky in RA, DEC .

Upvotes: 0

Hugh Bothwell
Hugh Bothwell

Reputation: 56714

For sake of reference, here is how I have arranged the octants:

enter image description here

from math import atan2, degrees

def toast_xy_to_latlon(x, y, size, inside=False):
    """
    Given TOAST x,y coordinates into a size*size HTM
    return lat, lon
      where -90 <= lat <= 90  (degrees North)
      and     0 <= lon <= 360 (degrees East)
    """
    half = size / 2

    # Figure out which quadrant (x, y) is in;
    #   save longitude offset and cast to first quadrant
    #   in new pixel-centered (u, v) coordinate system
    #   where 0.5 <= u,v <= half - 0.5 and (0,0) is North
    if x < half:
        if y < half:
            # 6,7
            lon_offset = 270
            u,v = half - x - 0.5, half - y - 0.5
        else:
            # 0,1
            lon_offset = 0
            u,v = y + 0.5 - half, half - x - 0.5
    else:
        if y < half:
            # 4,5
            lon_offset = 180
            u,v = half - y - 0.5, x + 0.5 - half
        else:
            # 2,3
            lon_offset = 90
            u, v = x + 0.5 - half, y + 0.5 - half

    # Find latitude
    lat = 90 * (1 - (u + v) / half)

    # Find longitude
    if u + v <= half:
        # Northern hemisphere
        lon = degrees(atan2(u, v))
    else:
        # Southern hemisphere
        lon = degrees(atan2(half - v, half - u))

    # Correct longitude offset -
    #   outside projection longitudes are measured CCW from left,
    #   inside projections CCW from right (rotated 180 degrees)
    if inside:
        lon_offset = (lon_offset + 180) % 360

    return lat, lon + lon_offset

Upvotes: 2

Related Questions