Plotting GPS track over Open Street Map background with Python

Question

I am trying to make one figure for each of my GPS hiking route where I plot the track on a map, very likely the one from open street map.

Extracting the data from the gpx file is not an issue and I have a pandas dataframe containing the coordinates (latitude and longitude), the time, and the altitude.

From the coordinates, I can download the tiles from open street map as described here: Easy OpenStreetMap tile displaying for Python

def plot_map(data1, zoom):

    lat_deg = data1['latitude'].min()
    delta_lat = data1['latitude'].max() - lat_deg
    lon_deg = data1['longitude'].min()
    delta_long = data1['longitude'].max() - lon_deg
    

    
    def deg2num(lat_deg, lon_deg, zoom):
      lat_rad = math.radians(lat_deg)
      n = 2.0 ** zoom
      xtile = int((lon_deg + 180.0) / 360.0 * n)
      ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
      return (xtile, ytile)
  
       
        
    def num2deg(xtile, ytile, zoom):
      n = 2.0 ** zoom
      lon_deg = xtile / n * 360.0 - 180.0
      lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
      lat_deg = math.degrees(lat_rad)
      return (lat_deg, lon_deg)
        
    
    

    headers = {"User-Agent":"Mozilla/5.0 (Macintosh; Intel Mac OS X 10_15_4) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/83.0.4103.97 Safari/537.36"}
    smurl = r"http://b.tile.openstreetmap.org/{0}/{1}/{2}.png"
    xmin, ymax =deg2num(lat_deg, lon_deg, zoom)
    xmax, ymin =deg2num(lat_deg + delta_lat, lon_deg + delta_long, zoom)
   
    Cluster = Image.new('RGB',((xmax-xmin+1)*256-1,(ymax-ymin+1)*256-1) ) 
    for xtile in range(xmin, xmax+1):
        for ytile in range(ymin,  ymax+1):
            try:
                imgurl = smurl.format(zoom, xtile, ytile)
                print("Opening: " + imgurl)
                imgstr = requests.get(imgurl, headers=headers)
                tile = Image.open(BytesIO(imgstr.content))
                Cluster.paste(tile, box = ((xtile-xmin)*256 ,  (ytile-ymin)*255))
            except: 
                print("Couldn't download image")
                tile = None

    fig = plt.figure()
    plt.imshow(Cluster)

However, the axis of this plot go from 0 to 256*n-1, with n the number of tiles in the x or y direction. Accordingly, I need to renormalize my coordinates from degrees into pixels.

Here, I am miserably failing, mostly because I cannot figure out the coordinates of the edges of my map.

Any suggestion on how to convert coordinates into pixels, given the number of tiles which depends on the zoom and on the track itself? Is there anything as simple as

import mplleaflet
mplleaflet.show()

but which saves the output to a jpg instead of opening an html file?

Thank you very much, Leonardo