将数字平铺到 lon./lat。在 EPSG:3395

Question

我需要将 tiles number 转换为 lon./lag in EPSG:3395 但我找不到解决方案。

我找到了 code for EPSG:4326 but i don't find a way to adapt it for EPSG:3395。

4326 代码（效果很好）：

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

我需要在 EPSG:3395 中将图块编号转换为 lon./lag，但我找不到可行的解决方案。我已经实现了一些基于 this answer.

的代码

当我尝试将 4326 度转换为 3395 度时：

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;

我收到以下错误消息：

    Unsupported operand type(s) for float and INT

更新：我已通过将 ^ 替换为 ** 来更正以下代码。

代码：

    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))

但是投影还是很奇怪。我认为问题出在这部分：

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

我不得不插入 math.radians，因为 math.tan 不喜欢度角。有什么想法吗？

Answer 1

您似乎想知道将坐标转换为 EPSG:3395 的公式。我不知道这是否是这个问题的正确位置，但这可能 help。

尝试乘幂时使用操作数 ** 而不是 ^。

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;

将数字平铺到 lon./lat。在 EPSG:3395

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

python

gis

tiles

projection