Tile numbers to lon./lat. in EPSG:3395

Question

I need to convert tiles number into lon./lag in EPSG:3395 but i don't find the solution.

I've found the code for EPSG:4326 but i don't find a way to adapt it for EPSG:3395.

Code for 4326 (it works well) :

$n = pow(2, $zoom);
$lon_deg = $xtile / $n * 360.0 - 180.0;
$lat_deg = rad2deg(atan(sinh(pi() * (1 - 2 * $ytile / $n))));

I need to convert tiles number into lon./lag in EPSG:3395 but i don't find a working solution. I have implemented some code based on this answer.

when I've tried to convert 4326 degrees into 3395 degrees:

def getVal(x, y, n):
    lon_deg = x / n * 360.0 - 180.0
    lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * y / n)))
    lat_deg = math.degrees(lat_rad)
    #Convert in 3395
    a = 6378137 #WGS84 semi-major axis
    b = 6356752.3142 #WGS84 semi-minor axis
    print(math.sqrt(1 - b^2 / a^2))
    e = math.sqrt(1 - b^2 / a^2) #ellipsoid eccentricity
    c = math.pow((1 - e*math.sin(latitude)) / (1 + e*math.sin(latitude)), e/2)
    lat_deg = a * ln(math.tan(math.pi/4 + lat_deg/2) * c)
    lon_deg = a * lon_deg; 

I obtain the following Error message:

    Unsupported operand type(s) for float and INT

Update: I have corrected to code to below by replacing ^ with **.

Code:

    def getVal(x, y, n):
        #Calcuate coordinates in 4326
        lon_deg = x / n * 360.0 - 180.0
        lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * y / n)))
        lat_deg = math.degrees(lat_rad)
        #Convert coordinates in 3395
        a = 6378137 #WGS84 semi-major axis
        b = 6356752.3142 #WGS84 semi-minor axis
        e = math.sqrt(1 - b**2 / a**2) #ellipsoid eccentricity
        c = math.pow((1 - e*math.sin(lat_deg)) / (1 + e*math.sin(lat_deg)), e/2)
        lon_deg = a * lon_deg;  
        lat_deg = a * math.log(math.tan(math.radians(math.pi/4 + lat_deg/2) * c))

But still the projection is strange. I assume the problem is on this part:

        lat_deg = a * math.log(math.tan(math.radians(math.pi/4 + lat_deg/2) 

I had to insert math.radians as math.tan doesn't like degree angle. Any idea?

Answer 1

Looks like you want to know the formulae to convert coordinates to EPSG:3395. I dont know if this is the right place for this question, but this might help.

When trying to raise to a power use the operand ** not ^.

def getVal(x, y, n):
    lon_deg = x / n * 360.0 - 180.0
    lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * y / n)))
    lat_deg = math.degrees(lat_rad)
    #Convert in 3395
    a = 6378137 #WGS84 semi-major axis
    b = 6356752.3142 #WGS84 semi-minor axis
    print(math.sqrt(1 - b**2 / a**2))
    e = math.sqrt(1 - b**2 / a**2) #ellipsoid eccentricity
    c = math.pow((1 - e*math.sin(latitude)) / (1 + e*math.sin(latitude)), e/2)
    lat_deg = a * ln(math.tan(math.pi/4 + lat_deg/2) * c)
    lon_deg = a * lon_deg;