Find pixel coordinate of world/geographic coordinate in tile

Question

I'm trying to use Mapbox Terrain RGB to get elevation for specific points in space. I used mercantile.tile to get the coordinates of the tile containing my point at zoom level 15, which for -43º, -22º (for simplicity sake) is 12454, 18527, then mercantile.xy to get the corresponding world coordinates: -4806237.7150042495, -2621281.2257876047.

Shouldn't the integer part of -4806237.7150042495 / 256 (tile size) equal the x coordinate of the tile containing the point, that is, 12454? If this calculation checked out I'd figure that I'm looking for the pixel column (x axis) corresponding to the decimal part of the result, like column 127(256 * 0,5) for 12454,5. However, the division results in -18774.366, (which is curiously close to the tile y coordinate, but it looks like a coincidence). What am I missing here?

As an alternative, I thought of using mercantile.bounds, assigning the first and last pixel columns to the westmost and eastmost longitudes, and finding my position with interpolation, but I wanted to check if I'm doing this the right/recommended way. I'm interested in point elevations, so everything said here goes for the Y axis as well.

Answer 1

Here's what I got so far:

def correct_altitude_mode(kml):
    with open(kml, "r+") as f:
        txt = f.read()
        if re.search("(?<=<altitudeMode>)relative(?=<\/altitudeMode>)", txt):
            lat = round(float(find_with_re("latitude", txt)), 5)
            lng = round(float(find_with_re("longitude", txt)), 5)
            alt = round(float(find_with_re("altitude", txt)), 5)
        z = 15
        tile = mercantile.tile(lng, lat, z)            
        westmost, southmost, eastmost, northmost = mercantile.bounds(tile)
        pixel_column = np.interp(lng, [westmost, eastmost], [0,256])
        pixel_row = np.interp(lat, [southmost, northmost], [256, 0])
        response = requests.get(f"https://api.mapbox.com/v4/mapbox.terrain-rgb/{z}/{tile.x}/{tile.y}.pngraw?access_token=pk.eyJ1IjoibWFydGltcGFzc29zIiwiYSI6ImNra3pmN2QxajBiYWUycW55N3E1dG1tcTEifQ.JFKSI85oP7M2gbeUTaUfQQ")
        buffer = BytesIO(response.content)
        tile_img = png.read_png_int(buffer)
        _,R,G,B = (tile_img[int(pixel_row), int(pixel_column)])
        print(tile_img[int(pixel_row), int(pixel_column)])
        height = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)
        print(f"R:{R},G:{G},B:{B}\n{height}")
        plt.hlines(pixel_row, 0.0, 256.0, colors="r")
        plt.vlines(pixel_column, 0.0, 256.0, colors="r")
        plt.imshow(tile_img)