How to figure out Latitude and Longitude from a x,y pixel on a Google Maps Static API image

Question

I am working on a function that is given the following parameters and returns latitude and longitude. The image is from the google maps static api. Google maps uses the Mercator projection. I am trying to figure out how to accurately get the latitude and longitude from the x,y pixel on the image. I ran multiple forms of calculations and they seem to be off slightly.

:param x, y: Pixel coordinates.
:param zoom: Zoom level of the image.
:param scale: Scale factor of the image.
:param image_center_lat_lon: Tuple of latitude and longitude at the image center.
:param image_size: Size of the image in pixels.

Here are two examples of functions that I have used that resulted in close but not accurate results

def pixel_to_lat_lon(x, y, zoom, scale, image_center_lat_lon, image_size):

    C = 640 * scale * (2 ** zoom)  # Total number of pixels in the map's width or height
    lat_center, lon_center = image_center_lat_lon
    image_size_scaled = (image_size[0] * scale, image_size[1] * scale)

    # Calculate pixel coordinates of the center of the image
    center_x = (lon_center + 180) / 360 * C
    center_y = (1 - math.log(math.tan(math.radians(lat_center)) + 1 / math.cos(math.radians(lat_center))) / math.pi) / 2 * C

    # Calculate the latitude and longitude from the pixel coordinates
    pixel_x = center_x - image_size_scaled[0] / 2 + x
    pixel_y = center_y - image_size_scaled[1] / 2 + y

    lon = pixel_x / C * 360 - 180
    lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * pixel_y / C)))
    lat = math.degrees(lat_rad)

    return lat, lon

Example 2:

def project(lat_lng, zoom, scale, tile_size):
    """Project latitude, longitude to world coordinates."""
    lat, lng = lat_lng
    siny = math.sin(lat * math.pi / 180.0)
    siny = min(max(siny, -0.9999), 0.9999)
    x = tile_size * (0.5 + lng / 360.0) * scale * (2 ** zoom)
    y = tile_size * (0.5 - math.log((1 + siny) / (1 - siny)) / (4 * math.pi)) * scale * (2 ** zoom)
    return x, y

def unproject(world_coordinate, zoom, scale, tile_size):
    """Unproject world coordinates to latitude, longitude."""
    x, y = world_coordinate
    lng = (x / (tile_size * scale * (2 ** zoom)) - 0.5) * 360.0
    n = math.pi - 2.0 * math.pi * y / (tile_size * scale * (2 ** zoom))
    lat = (180.0 / math.pi) * math.atan(0.5 * (math.exp(n) - math.exp(-n)))
    return lat, lng

def pixel_to_lat_lon(x, y, zoom, scale, image_center_lat_lon, image_size):
    tile_size = 640  
    center_x, center_y = project(image_center_lat_lon, zoom, scale, tile_size)
    
    # Calculate the center pixel of the image
    image_center_x = image_size[0] / 2
    image_center_y = image_size[1] / 2

    # Adjust x and y to be absolute pixel coordinates on the map
    abs_x = center_x + (x - image_center_x)
    abs_y = center_y + (y - image_center_y)
    
    # Convert absolute pixel coordinates back to latitude and longitude
    return unproject((abs_x, abs_y), zoom, scale, tile_size)

Here are some reference links to supporting documentation of the google maps api image:

Map and Tile Coordinates

Google Maps Static API Parameter Usage, helpful to understand zoom and scale

Test Case 1:

{
  "latitude": 25.64208985300285,
  "longitude": -80.37262545716628,
  "zoom": 18,
  "scale": 2,
  "size": 640,
  "pixelX": 767.4027099609375,
  "pixelY": 700.3953857421875
}

Test Case 2:

{
  "latitude": 25.63934327097157,
  "longitude": -80.35101258254728,
  "zoom": 18,
  "scale": 2,
  "size": 640,
  "pixelX": 247.67588806152344,
  "pixelY": 239.9282989501953
}

The results from the test case is a rough estimate. But should reflect the center of a pool according to the Google Maps UI, the latitude and longitude should be roughly around:

Test Case 1: 25.641942181931114, -80.37228858417026

Test Case 2: 25.640310082174018, -80.3520713140406

Answer 1

Solved using pyproj library with a Transformer to go between WGS84 and Google Maps Mercator

google_maps_proj = Proj('epsg:3857')
wgs84 = Proj('epsg:4326')

transformer = Transformer.from_proj(wgs84, google_maps_proj, always_xy=True)