How is it possible to calculate azimuth from a single lat lon cordinate with x, y offsets in meters

Question

Is it possible to calculate the azimuth from a single latitude and longitude and the destination's distance x and y from the mentioned coordinate?
Basically I have a single coordinate in WGS84, and I need to specify the next coordinate given the x, y distance in meters from this point.
I tried something like the following but it has a huge offset when distance increases:

def get_next_coordinate(initial_coordinate, dx_meters, dy_meters, is_up=False, is_right=True):
    lat, long = initial_coordinate
    # earth radius in km
    earth = 6378.137
    # 1 meter in degree
    coef = (1 / ((2 * math.pi / 360) * earth)) / 1000
    up_dir = 1 if is_up else -1
    right_dir = 1 if is_right else -1
    dx = (dx_meters * right_dir) * coef
    dy = (dy_meters * up_dir) * coef 
    new_lat = lat + dy
    new_long = long + dx / math.cos(lat * (math.pi/180))
    return new_lat, new_long

So, I thought I give geographiclib a try. it seems geograpghiclib.geodesicline.GeodesicLine has a SetDistance method that does exactly this; given a coordinate and a distance, it calculates the next coordinate. However it requires the azimuth for it to work.
Is it possible to somehow calculate the azimuth from what I have here? what are my options here?

More explanation:

To be more precise, I'm trying to calculate a latitude, longitude for any given pixel location in a map (which is an image). I have the latitude and longitude at map(0,0)
and I also have pixel scales, so I know how much space a pixel takes.
Using these two pieces of information, I can count pixels in x,y (height and width) and get the relative distances in meters.
My final step would be to calculate the coordinates at the given point.

Answer 1

After couple of days of searching and testing I finally managed to get this to work.
The first question was whether it was possible to calculate azimuth using an image's x,y pixel coordinate.
The answer is Yes.
Following the answer here, you can calculate the azimuth/bearing in a 2D image.
This is how the code looks in python:

def calculate_bearing(ax, ay, bx, by):
    """Computes the bearing in degrees from the point A(a1,a2) to
    the point B(b1,b2). Note that A and B are given in terms of
    screen coordinates.

    Args:
        ax (int): The x-coordinate of the first point defining a line.
        ay (int): The y-coordinate of the first point defining a line.
        bx (int): The x-coordinate of the second point defining a line.
        by (int): The y-coordinate of the second point defining a line.

    Returns:
        float: bearing in degrees
    """
    
    TWO_PI = math.pi * 2
    # if (a1 = b1 and a2 = b2) throw an error 
    theta = math.atan2(bx - ax, ay - by)
    if (theta < 0.0):
        theta += TWO_PI
    return math.degrees(theta) 

Also, there is this gist, where it calculates the direction of a line, in the same manner, suggested here.

# https://gist.github.com/aarondandy/604912
# for calculating the bearing/azimuth from image pixel coordinates
def calculate_line_direction (ax, ay, bx, by):
    """Calculates the north relative clockwise orientation of a line.

    Args:
        ax (int): The x-coordinate of the first point defining a line.
        ay (int): The y-coordinate of the first point defining a line.
        bx (int): The x-coordinate of the second point defining a line.
        by (int): The y-coordinate of the second point defining a line.

    Returns:
        float: The direction of the line in degree or nan on failure.
    """
    dx = bx - ax
    dy = by - ay
    m = math.sqrt((dx * dx) + (dy * dy))
    degree = math.degrees(math.atan2(dx / m, dy / m))
    return degree if m>0 else math.nan

Apart from this, I found out that calculating the bearing in my case can be done in a relatively simple way.

Since The azimuth is the heading measured clockwise from north, and our image is already aligned toward north, the point to the right will have 90 degree (clockwise) angle and the point below will have 180 degree angle.
We can hardcode these or calculate them like the method I just mentioned. this is a simple figure I came up with to show this.

    #  
    # |  -    
    # |      -
    # |         -
    # |           -
    # |             -
    # |               -
    # |                 -
    # |                  -
    # |                   -
    # |                    -
    # |                     -
    # point-a(0)--------point-b (90)
    # |                     - 
    # |                    -
    # |                   -
    # |                  -
    # |                -
    # |              -
    # |            -
    # |         -
    # |      -
    # |   -
    # point-c(180)

Before I continue, I want to talk about the pixel scale calculation. without a precise pixel scale calculation everything is in vein. Searching for how to calculate the pixel scale ends up in lots of results that dont work! like the one explained here :

# This doesn't work! don't use it!
const TILE_SIZE = 256;          // Tile size of base tile.
const CIRCUMFERENCE = 40075016; // Circumference at the equator in meters.
const scale = 1 << zoomLevel;   // Bitshifting is easy shorthand for Math.pow(2, zoomLevel)
const metersPerPixel = CIRCUMFERENCE / (TILE_SIZE * scale);

The thing is, if you use the snippet above, you'll get huge deviations based on where in the world you use it to calculate the pixel scale.
Basically when it comes to Mercator projection, it produces a different horizontal scale at different latitudes and thus pixel scale differs based on the latitudes.
What seems to get a very decent approximation is

# https://groups.google.com/g/google-maps-js-api-v3/c/hDRO4oHVSeM
# https://gis.stackexchange.com/a/127949/191141
def calculate_pixel_scale(lat, zoom):
    """Calculates and returns the pixel scale in meters at the specified zoom level.
        This function shows what a pixel represents in real world. using this you can get
        the real world estimate of distance traversed, or object sizes, etc in meters. 
        This is based on the GoogleMap's specifications (zoom, tile size, etc). 

    Args:
        lat (float): the latitude at which you want to calculate the pixel scale
        zoom (int): the zoom level at which you want to calculate the pixel scale

    Returns:
        float: pixel scale
    """
    scale = 156543.03392 * math.cos(lat * math.pi / 180) / math.pow(2, zoom)
    return scale

Now, somehow in order to get this, and calculate the next latitude/longitude, I need to have the latitude of that location! it seems like a circle now!
In order to break this cycle, we can use google tile coordinates. we know that at different zooms, we have different number of tiles that constitute our map.
These tiles are all 256x256 in size, so if we know a single tile at our specified zoom level in our map, we can calculate any other tile coordinates in our map (as they are 256x256, its a simple math)
So if we can get the tile coordinate at our map(0,0) pixel coordinate, we be fine. there are many ways to do this, but one simple way would be to calculate tile coordinate from a given latitude, longitude. This can be done like this:

# calculate tile coordinate based on lat/lon
# https://github.com/Coderx7/google-map-downloader/blob/master/downloader_1.1.py#L110
# note read the note!
def calculate_tile_coordinate(latitude, longitude, zoom):
    """Get google-style tile cooridinate from geographical coordinate
    Note: 
    Note that this will return a tile coordinate. and a tile coordinate
    is located at (0,0) or origin of a tile. 
    all tiles are 256x256. there are as many gps locations inside of a tile 
    that when used with this function will return the same exact tile coordinate
    (as they all reside in that tile obviously) in order to get the displacement 
    of your gps location, after calculating the tile coordinate, calculate tile's 
    (0,0) gps location from the newly achieved tile coordinate. then you have gps location
    at (0,0) and now can calculate howfar your initial gps location is from the 
    origin of the tile. 

    Args:
        latitude (float): Latitude
        longitude (float): Longittude
        zoom (int): zoom

    Raises:
        TypeError: [description]
        TypeError: [description]

    Returns:
        tuple(int,int): returns tile coordinate in the form of (y,x)
    """
    isnum = lambda x: isinstance(x, int) or isinstance(x, float)
    if not (isnum(longitude) and isnum(latitude)):
        raise TypeError("latitude and longitude must be int or float!")

    if not isinstance(zoom, int) or zoom < 0 or zoom > 22:
        raise TypeError("zoom must be int and between 0 to 22.")

    if longitude < 0:
        longitude = 180 + longitude
    else:
        longitude += 180
    longitude /= 360  # make longitude to (0,1)

    latitude = 85.0511287798 if latitude > 85.0511287798 else latitude
    latitude = -85.0511287798 if latitude < -85.0511287798 else latitude
    latitude = math.log(math.tan((90 + latitude) * math.pi / 360)) / (math.pi / 180)
    latitude /= 180  # make latitude to (-1,1)
    latitude = 1 - (latitude + 1) / 2  # make latitude to (0,1) and left top is 0-point

    num = 2 ** zoom
    y = math.floor(latitude * num)
    x = math.floor(longitude * num)
    return y,x

Now that we have the tile coordinate, we can check at which tile our desired pixel resides, then we get that tile's gps coordinate, now we can use it's latitude to calculate the pixel scale, and then calculate the offset inside that tile, to get the exact gps location of specified pixel coordinate.
also now that we can have any tile number, by using the tile to the right and the tile below the target tile, we can calculate exact pixel scale for x, and y axis respectively(shown below).

This is my quick and dirty code that I have been playing wit,but before that, for better visualization this might help:

def get_distance(gps_coordinate1, gps_coordinate2):
    """Calculates the distance between to gps coordinates in meters.

    Args:
        coord1 (tuple(float, float)): latitude and longitude of point A
        coord2 (tuple(float, float)): latitude and longitude of point B

    Returns:
        float: distance is meters
    """
    return Geodesic.WGS84.Inverse(gps_coordinate1[0],
                                  gps_coordinate1[1], 
                                  gps_coordinate2[0], 
                                  gps_coordinate2[1])['s12']

def calculate_latitue_longitude_from_tile(row, column, zoom_level=18):
    """ 
        Given a tile coordinate, calculates the latitude and longitude at 
        the specified zoom level.

    Args:
        row (int): tile's row number. 
        column (int): tile's column number. 
        zoom_level (int, optional): the zoom at which level, map is collected. Defaults to 18.

    Returns:
        [type]: [description]
    """
    lon = column / math.pow(2.0, zoom_level) * 360.0 - 180
    n = math.pi - (2.0 * math.pi * row) / math.pow(2.0, zoom_level)
    lat = -((math.atan(math.sinh(n))) * -1) / math.pi * 180
    return lat, lon

def calculate_gps_from_pixel_offsets(x1=256, y1=0, offset_x=256, offset_y=0, target_tile=None, use_offsets_for_bearing=False, zoom=18):
    """
        calculates a new coordinate given an initial coordinate and pixel values x,y 
        In order to calculate a new coordinate from an existing one with a given distance 
        in meters, we need to calculate the bearing/azimuth. that is the direction at which
        we are headed needs to be known. for this we can calculate the angle between x,y 
        in our image. 

        the image we just talked about is our map. The map that was 
        from google earth. Knowing the tile numbers, we can calculate that locations coordinate
        since tiles are 256x256, and we created our map from these tiles, the (0,0) of our map
        corresponds to a tile number, we convert this tile coordinates to gps coorinate
        and we have our first coordinate ready. now we need to specify at what pixel location 
        we want to calculate the respective gps location.

        Suppose at (512,768), we calculate the angle between these pixels (row,col) with respec to
        north (googlemap is aligned to north) since we have coordinate at 0,0 of the map, other 
        locations will be on right(east) and down(south) direction. so calculating the azimuth from 
        something like (0,256) and (256,0) would apply to all other values of x,y in our map. 
        
        So we simply set x1,andy1 to 256 and 0 respectively. (0,256) will give us the azimuth for 
        x axis, and (256,0) will give us the azimuth for y axis, and using these two azimuths we 
        calculate two new locations/coordinates we use latitude from the one that uses y-azimuth 
        and use the longitude that was calculated using x-azimuth. this will give us the final 
        coordinate that we want. 

        as you can see, we also need the distance, offset_x and offset_y, they represent the 
        distance from our target tile pixel location. since they are in pixels
        we need to calculate the pixel scale at the zoom level we took the image, this will allow us
        to multiply the pixel offsets with this scale and get the metric distance. we can then
        use this metric distance to reach from the initial coordinate to the final one that we are after

    Args:
        x1 (int, optional): target x-pixel coord, we want to get its gps coord. Defaults to 256.
        y1 (int, optional): target y-pixel coord, we want to get its gps coord. Defaults to 0.
        offset_x (int, optional): pixel distance in x axis from tile's origin. Defaults to 256. it can be a value from 0-256
        offset_y (int, optional): pixel distance in y axis from tile's origin. Defaults to 0. it can be a value from 0-256
        target_tile (tuple(float,float), optional): the target tile coordinate that's used for calculating other coordinates for each pixel. Defaults to None.
        use_offsets_for_bearing (bool, optional): whether to use pixel coordinates for calculating azimuth instead of x1,y1. Defaults to False.
    """
    
    if use_offsets_for_bearing:
        x1, y1 = offset_x, offset_y

    # calculate azimuth/bearing for x,y individually
    azimuth_x = calculate_bearing(0, 0, x1, 0)
    azimuth_y = calculate_bearing(0, 0, 0, y1)

    # calculate the gps coordinate
    target_tile_gps = calculate_latitue_longitude_from_tile(target_tile[0], target_tile[1], zoom)

    # pixel scale at zoom 18
    # This calculates the azimuth, and is fairly accurate, however
    # the better way for calculating the azimuth is by using a separate
    # scaler for x,y respectively, this can be achieved by calculating the 
    # azimuth from 3 gps locations, one to the right and one down. 
    # these gps locations can be easily calculate by having one tile which
    # is our target tile, the tile that contains our pixel coordinate x,y
    # we then can simply calculate the other two tiles using simple math 
    # simply (tile_row, tile_col+1) and (tile_row+1, tile_col) will give us
    # the tiles, which we can then calculate the gps using google's formula
    # given before. this is done in our second method
    pixel_scale = calculate_pixel_scale(target_tile_gps[0], 18)

    # calculating the offset in meters based on pixel scale
    dsx_meters = offset_x * pixel_scale
    dsy_meters = offset_y * pixel_scale
    
    # if we set x1 and y1 to actual pixel locations(dsx,dsy), we need to remove the 180
    if not use_offsets_for_bearing:
        azimuth_x = 180 - azimuth_x 
        azimuth_y = 180 - azimuth_y 
   
    # to right/east direction
    g = Geodesic.WGS84.Direct(target_tile_gps[0], target_tile_gps[1], azimuth_x, dsx_meters)
    # to down/south direction
    g2 = Geodesic.WGS84.Direct(target_tile_gps[0], target_tile_gps[1], azimuth_y, dsy_meters)
    
    final_gps = g2['lat2'], g['lon2']

    print(f'\n---------------Method I-------------')
    print(f'distance offset in x axis: {dsx_meters:.2f}m.\ndistance offset in y axis: {dsy_meters:.2f}m.')
    print(f"GPS based on dsx: ({g['lat2']:.8f}, {g['lon2']:.8f}) az1: {azimuth_x:.2f}\N{DEGREE SIGN}")
    print(f"GPS based on dsy: ({g2['lat2']:.8f}, {g2['lon2']:.8f}) az2: {azimuth_y:.2f}\N{DEGREE SIGN}")
    print(f"GPS based on dsx-dsy: ({g2['lat2']:.8f}, {g['lon2']:.8f}).")
    print(f'Tile GPS: {target_tile_gps}')
    print(f"Final GPS:{final_gps}")
    print(f'Distance from tile gps to calculared gps location: {get_distance(target_tile_gps, final_gps):.2f}m.')

    return final_gps

def calculate_gps_from_pixel_offsets_accurate(target_tile, offset_x, offset_y, zoom=18):
   
    # calculate new tile's coordinate
    row, column = target_tile
    tile_right = (row, column+1)
    tile_down = (row+1, column)

    # calculate gps coordinates
    tile_gps = calculate_latitue_longitude_from_tile(row, column, zoom)
    tile_right_gps = calculate_latitue_longitude_from_tile(row, column+1, zoom)
    tile_down_gps = calculate_latitue_longitude_from_tile(row+1, column, zoom)

    # Define the path from 1 to 2
    lx = Geodesic.WGS84.InverseLine(tile_gps[0], tile_gps[1], tile_right_gps[0], tile_right_gps[1])
    ly = Geodesic.WGS84.InverseLine(tile_gps[0], tile_gps[1], tile_down_gps[0], tile_down_gps[1])
    
    # calculating pixel scale:
    # between tiles, there is a 256x256 distance in (x,y)
    pixel_scale_x = lx.s13/256
    pixel_scale_y = ly.s13/256

    # distance in meters in x axis
    dx_meters = offset_x*pixel_scale_x
    # distance in meters in y axis
    dy_meters = offset_y*pixel_scale_y
    
    # https://gis.stackexchange.com/questions/349458/
    direction_x = Geodesic.WGS84.Direct(tile_gps[0], tile_gps[1], lx.azi1, dx_meters)
    direction_y = Geodesic.WGS84.Direct(tile_gps[0], tile_gps[1], ly.azi1, dy_meters)
    
    final_gps = direction_y['lat2'], direction_x['lon2']

    print(f'\n---------------Method II-------------')
    print(f"Final GPS: {final_gps}")
    print(f'azi1_x: {lx.azi1}')
    print(f'azi1_y: {ly.azi1}')
    print(f'pixel-scale-x: {pixel_scale_x}')
    print(f'pixel-scale-y: {pixel_scale_y}')
    print(f'dx in meters: {dx_meters:.2f}m')
    print(f'dy in meters: {dy_meters:.2f}m')

    return final_gps

def find_tile_using_pixel_xy(first_tile, pixel_coordxy, zoom=18):
    """
        Given a pixel coordinate, calculate the target tile and (dx, dy) for tile-level coordinate.   
        dx, dy are tile level coordinates, which are in the range 0-256 representing, x, y offset `inside` of the given tile.

    Args:
        first_tile (tuple(int,int)): The first tile coordinate in our map image (0,0)(the upper right most corner of our image)
        pixel_coord (tuple(in,int)): The pixel coordinate for which we want the gps (with respect to initial tile location)
        zoom (int): the zoom level at which the tiles are specified (Default = 18)
    """

    row_init, col_init = first_tile
    x,y = pixel_coordxy

    # Calculate the steps needed to get from initial tile
    # to get to the target tile where our pixels reside in
    # since tiles are 256x256, we simply divide by 256 and this
    # will give us the steps needed for rows and columns 
    # example: x = 1022, y = 769 
    step_x = int(x/256)
    step_y = int(y/256)
    
    # Since we need to know where these global pixel coordinate reside exactly 
    # in our 256by256 tile internal coordinate. 
    # we do this by calculating the offset by taking the remainder from the division by 256 
    dx = abs(x % 256)
    dy = abs(y % 256) 

    selected_tile_row = row_init + step_y
    selected_tile_col = col_init + step_x
    selected_tile = (selected_tile_row, selected_tile_col)

    print(f'Pixel coordinate query: {pixel_coordxy}')
    print(f"Relative offsets inside tile: dx:{int(dx)} dy:{int(dy)}")    
    print(f'First tile:{first_tile} xy-step:({step_x},{step_y})')
    print(f'Selected tile:{selected_tile}')

    return selected_tile, (dx, dy)

def find_gps_using_pixel_xy(first_tile, pixel_coordxy, zoom=18):

    selected_tile, (dx, dy) = find_tile_using_pixel_xy(first_tile=first_tile, 
                                                       pixel_coordxy=pixel_coordxy,
                                                       zoom=zoom)
    # calculate the gps coordinate
    selected_tile_gps = calculate_latitue_longitude_from_tile(selected_tile[0], selected_tile[1], zoom)

    # first method - using approximate scale for x,y
    final_gps1 = calculate_gps_from_pixel_offsets(target_tile=selected_tile, 
                                                  offset_x=dx, 
                                                  offset_y=dy,
                                                  use_offsets_for_bearing=False,
                                                  zoom = zoom)
    
    # second method(the most accurate) - using separate pixel scales for x,y
    final_gps2 = calculate_gps_from_pixel_offsets_accurate(target_tile=selected_tile, 
                                                           offset_x=dx, 
                                                           offset_y=dy)
    
    print(f'Selected GPS:      {selected_tile_gps}')
    print(f'Final GPS-I:       {final_gps1}')
    print(f'Final GPS-II(acc): {final_gps2}')
    print(f'Distance from selected tile(GPS-I):  {get_distance(selected_tile_gps, final_gps1):.2f}m.')    
    print(f'Distance from selected tile(GPS-II): {get_distance(selected_tile_gps, final_gps2):.2f}m.')    
    print(f'-----------------------------------')

    return final_gps1, final_gps2

test:

init_gps = (38.213367074438175, -98.20816040039062)  
# when caclulating the tile coordinate using gps, make sure you account 
# for the displacement. see my note for this function above
tile = calculate_tile_coordinate(init_gps[0], init_gps[1], 18) # (100919, 59559)
tile_gps = calculate_latitue_longitude_from_tile(tile[0], tile[1], 18)
# they should match and indeed they do!
print(f'tile {tile} gps:{init_gps}')
print(f'tile {tile} gps:{tile_gps}')
# now lets test
# get the gps location at pixel coordinate (50,100)
# remember this pixel coordinate is relative to our first_tile
# we assumed that pixel(0,0) of our image is where our x0,y0
# of our first tile is.
find_gps_using_pixel_xy(tile, (50,100), 18)

output:

tile (100919, 59559) gps:(38.213367074438175, -98.20816040039062)
tile (100919, 59559) gps:(38.213367074438175, -98.20816040039062)
Pixel coordinate query: (50, 100)
Relative offsets inside tile: dx:50 dy:100
First tile:(100919, 59559) xy-step:(0,0)
Selected tile:(100919, 59559)

---------------Method I-------------
distance offset in x axis: 23.46m.
distance offset in y axis: 46.92m.
GPS based on dsx: (38.21336707, -98.20789252) az1: 90.00°
GPS based on dsy: (38.21294437, -98.20816040) az2: 180.00°
GPS based on dsx-dsy: (38.21294437, -98.20789252).
Tile GPS: (38.213367074438175, -98.20816040039062)
Final GPS:(38.212944374143305, -98.20789252325251)
Distance from tile gps to calculared gps location: 52.46m.

---------------Method II-------------
Final GPS: (38.21294558225481, -98.20789217948914)
azi1_x: 89.9995752467753
azi1_y: 180.0
pixel-scale-x: 0.46980161807460397
pixel-scale-y: 0.46785849030688803
dx in meters: 23.49m
dy in meters: 46.79m
--------------------------------------
Target Tile GPS:   (38.213367074438175, -98.20816040039062)
Final GPS-I:       (38.212944374143305, -98.20789252325251)
Final GPS-II(acc): (38.21294558225481, -98.20789217948914)
Distance from selected tile(GPS-I):  52.46m.
Distance from selected tile(GPS-II): 52.35m.