Why is my OSM map rendering squished in Love2D?

I'm trying to render an OpenStreetMap (.osm) file in Love2D without using external libraries. The issue is that the displayed map appears vertically squished instead of maintaining its correct proportions.

Here is how I process the .osm file:

1. parse the XML to extract nodes and ways.
2. convert latitude and longitude directly to (x, y) coordinates.
3. compute a scale factor based on the bounding box of the map and the screen dimensions (same for X and Y axis)
4. draw the map using love.graphics.line() as precalculated line.

Problem: The map appears compressed vertically instead of maintaining the correct aspect ratio.

The code:

local filenames = {
    'map-55.osm', --  -- 54.92349, -2.96367
    'map-0.osm', --  -- -0.697995, 10.243917
}

-- osm
local osm = {}

function osm.parseOSM(xml)
    local nodes = {}
    local nodesHash = {}

    local nodeIndex = 0
    local ways = {}
    for wayStr in xml:gmatch('<way.->.-</way>') do
        local nodeIndices = {}
        for id in wayStr:gmatch('<nd ref="(.-)"') do -- id or nd
            if nodesHash[id] then
                table.insert (nodeIndices, nodesHash[id])
            else
                nodeIndex = nodeIndex + 1
                nodesHash[id] = nodeIndex
                table.insert (nodeIndices, nodeIndex)
            end
        end
--      print ('#nodeIndices', #nodeIndices)
        local way = {nodeIndices = nodeIndices, line = {}}
        table.insert (ways, way)
    end

    for nodeStr in xml:gmatch('<node.-/>') do
        local id = nodeStr:match('id="(.-)"')
        local x = tonumber(nodeStr:match('lon="(.-)"')) -- X
        local y = tonumber(nodeStr:match('lat="(.-)"')) -- Y
        if nodesHash[id] then
            nodes[nodesHash[id]] = {x=x, y=-y}
        else

        end
    end

    local minLat, minLon, maxLat, maxLon = xml:match('<bounds minlat="(.-)" minlon="(.-)" maxlat="(.-)" maxlon="(.-)"/>')

    local minX = tonumber(minLon) -- x
    local maxX = tonumber(maxLon) -- x
    local minY = tonumber(minLat) -- y
    local maxY = tonumber(maxLat) -- y

    local dx = maxX-minX
    local dy = maxY-minY

    local bounds = {
        minX=minX, 
        maxX=maxX, 
        dx=dx,
        midX=minX + dx/2, 
        minY=minY, 
        maxY=maxY, 
        dy=dy,
        midY=minY + dy/2, 
    }

    return {nodes = nodes, ways = ways, bounds = bounds}
end

function osm.updateScale(mapData)
    local bounds = mapData.bounds
    local screenWidth, screenHeight = love.graphics.getDimensions ()
    mapData.scale = math.min(screenWidth / bounds.dx, screenHeight / bounds.dy)
    print ('mapData.scale', mapData.scale)

    local ways = mapData.ways
    local nodes = mapData.nodes
    for i, way in ipairs (ways) do
        way.line = {}
        for j, nodeIndex in ipairs (way.nodeIndices) do
            local node = nodes[nodeIndex]
            local x = node.x * mapData.scale
            local y = node.y * mapData.scale
            print ('x: '..x, 'y: '..y)
            table.insert (way.line, x)
            table.insert (way.line, y)
        end
    end
end

local mapDataSet = {}

function love.load()
    for i, filename in ipairs (filenames) do
        local file = love.filesystem.read(filename)
        if file then
            local mapData = osm.parseOSM(file)
            osm.updateScale(mapData)
            mapDataSet[i] = mapData
        end
    end
    mapData = mapDataSet[1]
end

function love.draw()
    local tx = mapData.bounds.midX*mapData.scale - 400
    local ty = mapData.bounds.midY*mapData.scale + 300
    love.graphics.translate (-tx, ty)
    for _, way in pairs(mapData.ways) do
        love.graphics.line(way.line)
    end
end

At higher latitudes (lat 55°), the geometry looks compressed vertically, as if the map is being flattened:

Near the equator (lat ≈ 0°), the shapes appear correct:

How can I correctly transform latitude values to maintain the correct proportions across the whole map?

Edited: I've used the Mercator projection:

local function mercatorY(lat)
    local radLat = math.rad(lat)
    local merY = math.deg(math.atan(math.sinh(radLat)))
    return merY
end

but it was not enough, I need extra curved factor with magic value to make the 55° round again:

local function mercatorY(lat)
    local radLat = math.rad(lat)
    local merY = math.deg(math.atan(math.sinh(radLat)))
    local stretchY = 1 + (1 - math.cos(radLat)) * 1.1
    return merY * stretchY
end

Than I've tried to make check the other lat (lat = 40°) and it was ellipse too, but too high:

    'map-40.osm', --  -- 39.875064, 20.027293