Unity 3d - Web Mercator image to Sphere Shader

I'm using a Web Mercator image to cover a sphere. My shader takes a plane with an image and turns it into sphere. The only issue is that The resulting sphere ends up with countries stretched (like the united states).

I've figured out that I can use an equlateral image of earth to get the desired effect of non-stretched countries

Question

For my project I only have web mercator imagery and I've been struggling with the math for getting my shader to show countries at their correct scale. How can I transform mercator lat lon to equilateral lat lon for writing to my shader ?

NOTE

Everything I would need seems to be on this question about mercator projection to equirectangular but for whatever reason it's just not clicking.

Some Code

plane script

using System.Collections;
using System.Collections.Generic;
using UnityEngine;

public class SquareBender : MonoBehaviour
{
    // Start is called before the first frame update
    void Start()
    {
        Vector3Int tileIndex = new Vector3Int(0, 0, 0);
        Mesh mesh = GetComponent<MeshFilter>().mesh;
        this.SetUpTileLonLats(mesh, tileIndex);

       GetComponent<Renderer>().material.SetFloat("SphereRadius", 50);
    }

    // tileIndex is column/row/zoom of current tile
    // uv is relative postion within tile
    //   (0,0) for bottom left, (1,1) top right
    Vector2 GetLonLatOfVertex(Vector3Int tileIndex, Vector2 uv)
    {

        float lon = uv.x * 360 - 180;
        // float lat = uv.y * 180 - 90;
        float lat = uv.y * 168 - 84;

        float lonRad = lon / 180 * Mathf.PI;//uv.x * Mathf.PI * 2 - Mathf.PI; 
        float latRad = lat / 180 * Mathf.PI;//uv.y * Mathf.PI - Mathf.PI / 2;

        float theta = lonRad;
        float phi = Mathf.Log(Mathf.Tan(Mathf.PI/4 + latRad/2));

        Debug.Log($"{uv.x} {uv.y} -- {lon} {lat} -- {lonRad} {latRad} -- {theta} {phi}");

        // Use tileIndex and uv to calculate lon, lat (in RADIANS)
        // Exactly how you could do this depends on your tiling API...
    
        return new Vector2(theta, phi);
    }

    // Call after plane mesh is created, and any additional vertices/uvs are set
    // tileIndex is column/row/zoom of current tile
    void SetUpTileLonLats(Mesh mesh, Vector3Int tileIndex)
    {
        Vector2[] uvs = mesh.uv;
        Vector2[] lonLats= new Vector2[uvs.Length];
        
        for (int i = 0; i < lonLats.Length; i++)
        {
            lonLats[i] = GetLonLatOfVertex(tileIndex, uvs[i]);
        }
        
        mesh.uv2 = lonLats;
    }

}

shader

Shader "Custom/SquareBender" {
    Properties{
        _MainTex("Tex", 2D) = "" {}
        _SphereCenter("SphereCenter", Vector) = (0, 0, 0, 1)
        _SphereRadius("SphereRadius", Float) = 50
    }

    SubShader{
        Cull off // for doublesized texture @jkr todo: disable for prod
        Pass {
            CGPROGRAM
                #pragma vertex vert
                #pragma fragment frag

                #include "UnityCG.cginc"

                struct appdata {
                   float2 uv     : TEXCOORD0;
                   float2 lonLat : TEXCOORD1;
                };

                struct v2f
                {
                    float4 pos  : SV_POSITION;
                    float3 norm : NORMAL;
                    float2 uv   : TEXCOORD0;
                };

                float4 _SphereCenter;
                float _SphereRadius;

                v2f vert(appdata v)
                {
                    v2f o;
                    float lon = v.lonLat.x;
                    float lat = v.lonLat.y;

                    _SphereRadius = 40;

                    fixed4 posOffsetWorld = fixed4(
                        _SphereRadius*cos(lat)*cos(lon),
                        _SphereRadius*sin(lat),
                        _SphereRadius*cos(lat)*sin(lon), 0);

                    float4 posObj = mul(unity_WorldToObject,
                            posOffsetWorld + _SphereCenter);

                    o.pos = UnityObjectToClipPos(posObj);
                    o.uv = v.uv;
                    o.norm = mul(unity_WorldToObject, posOffsetWorld);
                    return o;
                }

                sampler2D _MainTex;

                float4 frag(v2f IN) : COLOR
                {
                    fixed4 col = tex2D(_MainTex, IN.uv);
                    return col;
                }
            ENDCG
        }
    }
    FallBack "VertexLit"
}

** EDIT **

Shout out to @Ruzihm for his shader contribution from this answer about wrapping map tiles around a sphere

Solution

This question is related to a much larger tile-based earth question which I still haven't solved (at the time of writing, has since been resolved)

BUT

I was able to figure out how to solve this sub-question by using math from a potentially helpful answer I mentioned earlier in the OP

The Solution

I took the some of the shader code from @Pluto's answer and merged it in with my current shader. I assigned a web mercator image to a plane that also had this shader attached to it. The default "Projection" shader param is 0 so everything is already set to convert the mercator image to equirectangular and viola~ the image is rendered as equirectangular on a sphere.

MercatorBender.shader

Shader "Custom/MercatorBender" {
    Properties{
        _MainTex("Tex", 2D) = "" {}
        _SphereCenter("SphereCenter", Vector) = (0, 0, 0, 1)
        _SphereRadius("SphereRadius", Float) = 5

        [Enum(Equirectangular,0,Azimuthal,1)]
        _Azimuthal("Projection", float) = 0
    }

    SubShader{
        Cull off // for doublesized texture @jkr todo: disable for prod
        Pass {
            CGPROGRAM
                #pragma vertex vert
                #pragma fragment frag

                #include "UnityCG.cginc"

                struct appdata {
                   float2 uv     : TEXCOORD0;
                   float2 lonLat : TEXCOORD1;
                };

                struct v2f
                {
                    float4 pos  : SV_POSITION;
                    float3 norm : NORMAL;
                    float2 uv   : TEXCOORD0;
                };

                float4 _SphereCenter;
                float _SphereRadius;

                v2f vert(appdata v)
                {
                    v2f o;
                    float lon = v.lonLat.x;
                    float lat = v.lonLat.y;

                    _SphereRadius = 40;

                    fixed4 posOffsetWorld = fixed4(
                        _SphereRadius*cos(lat)*cos(lon),
                        _SphereRadius*sin(lat),
                        _SphereRadius*cos(lat)*sin(lon), 0);

                    float4 posObj = mul(unity_WorldToObject,
                            posOffsetWorld + _SphereCenter);

                    o.pos = UnityObjectToClipPos(posObj);
                    o.uv = v.uv;
                    o.norm = mul(unity_WorldToObject, posOffsetWorld);
                    return o;
                }

                sampler2D _MainTex;
                float _Azimuthal;

                // float4 frag(v2f IN) : COLOR
                // {
                //     fixed4 col = tex2D(_MainTex, IN.uv);
                //     return col;
                // }

#define PI 3.141592653589793238462f
#define PI2 6.283185307179586476924f

                float2 uvToEquirectangular(float2 uv) {
                    float lat = (uv.x) * PI2;   // from 0 to 2PI
                    float lon = (uv.y - .5f) * PI;  // from -PI to PI
                    return float2(lat, lon);
                }

                float2 uvAsAzimuthalToEquirectangular(float2 uv) {                  
                    float2 coord = (uv - .5) * 4; 

                    float radius = length(coord);
                    float angle = atan2(coord.y, coord.x) + PI;

                    //formula from https://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection
                    float lat = angle;
                    float lon = 2 * acos(radius / 2.) - PI / 2;
                    return float2(lat, lon);
                }       

                fixed4 frag(v2f i) : SV_Target
                {
                    // get equirectangular coordinates
                    float2 coord = _Azimuthal ? uvAsAzimuthalToEquirectangular(i.uv) : uvToEquirectangular(i.uv);

                    // equirectangular to mercator
                    float x = coord.x;
                    float y = log(tan(PI / 4. + coord.y / 2.));
                    // brin x,y into [0,1] range
                    x = x / PI2;
                    y = (y + PI) / PI2;                 

                    fixed4 col = tex2D(_MainTex, float2(x,y));

                    // just to make it look nicer
                    col = _Azimuthal && length(i.uv*2-1) > 1 ? 1 : col;

                    return col;
                }
            ENDCG
        }
    }
    FallBack "VertexLit"
}