如何在 python 或 Javascript 中从 UTM 转换为 LatLng

IT技术 javascript python gis arcgis-js-api proj4js
2021-02-27 15:47:20

我有一堆带有 UTM 形式坐标的文件。对于每个坐标,我都有东距、北距和区域。我需要将其转换为 LatLng 以与 Google Map API 一起使用以在地图中显示信息。

我发现一些在线计算器可以做到这一点,但没有实际的代码或库。http://trac.osgeo.org/proj4js/是 Javascript 的投影库,但查看演示它不包括 UTM 投影。

我对整个 GIS 领域还很陌生,所以我想要的是一些东西:

(lat,lng) = transform(easting, northing, zone)
6个回答

我最终从 IBM 找到了解决它的 Java 代码:http : //www.ibm.com/developerworks/java/library/j-coordconvert/index.html

仅供参考,这是我需要的方法的python实现:

import math

def utmToLatLng(zone, easting, northing, northernHemisphere=True):
    if not northernHemisphere:
        northing = 10000000 - northing

    a = 6378137
    e = 0.081819191
    e1sq = 0.006739497
    k0 = 0.9996

    arc = northing / k0
    mu = arc / (a * (1 - math.pow(e, 2) / 4.0 - 3 * math.pow(e, 4) / 64.0 - 5 * math.pow(e, 6) / 256.0))

    ei = (1 - math.pow((1 - e * e), (1 / 2.0))) / (1 + math.pow((1 - e * e), (1 / 2.0)))

    ca = 3 * ei / 2 - 27 * math.pow(ei, 3) / 32.0

    cb = 21 * math.pow(ei, 2) / 16 - 55 * math.pow(ei, 4) / 32
    cc = 151 * math.pow(ei, 3) / 96
    cd = 1097 * math.pow(ei, 4) / 512
    phi1 = mu + ca * math.sin(2 * mu) + cb * math.sin(4 * mu) + cc * math.sin(6 * mu) + cd * math.sin(8 * mu)

    n0 = a / math.pow((1 - math.pow((e * math.sin(phi1)), 2)), (1 / 2.0))

    r0 = a * (1 - e * e) / math.pow((1 - math.pow((e * math.sin(phi1)), 2)), (3 / 2.0))
    fact1 = n0 * math.tan(phi1) / r0

    _a1 = 500000 - easting
    dd0 = _a1 / (n0 * k0)
    fact2 = dd0 * dd0 / 2

    t0 = math.pow(math.tan(phi1), 2)
    Q0 = e1sq * math.pow(math.cos(phi1), 2)
    fact3 = (5 + 3 * t0 + 10 * Q0 - 4 * Q0 * Q0 - 9 * e1sq) * math.pow(dd0, 4) / 24

    fact4 = (61 + 90 * t0 + 298 * Q0 + 45 * t0 * t0 - 252 * e1sq - 3 * Q0 * Q0) * math.pow(dd0, 6) / 720

    lof1 = _a1 / (n0 * k0)
    lof2 = (1 + 2 * t0 + Q0) * math.pow(dd0, 3) / 6.0
    lof3 = (5 - 2 * Q0 + 28 * t0 - 3 * math.pow(Q0, 2) + 8 * e1sq + 24 * math.pow(t0, 2)) * math.pow(dd0, 5) / 120
    _a2 = (lof1 - lof2 + lof3) / math.cos(phi1)
    _a3 = _a2 * 180 / math.pi

    latitude = 180 * (phi1 - fact1 * (fact2 + fact3 + fact4)) / math.pi

    if not northernHemisphere:
        latitude = -latitude

    longitude = ((zone > 0) and (6 * zone - 183.0) or 3.0) - _a3

    return (latitude, longitude)

在这里,我认为这是一些简单的easting*x+zone*y事情。

@Franki:因为 UTM 坐标是相对于预定义区域的。维基百科上的 UTM
2021-04-26 15:47:20
我不怪你认为这很简单。我上周从经纬度转换为米勒投影的东/北,我还在里面庆祝。
2021-05-10 15:47:20
你救了我的一天,我为这样的事情找了将近 2 天,而且效果很好,我做了你脚本的 javascript 版本,如果其他人觉得它有帮助,我会把它作为答案......非常感谢你
2021-05-10 15:47:20
为什么您的功能需要区域?你不应该只需要东移,北移,半球吗?
2021-05-18 15:47:20

我找到的是以下网站:http : //home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html 它有一个 javascript 转换器,你应该检查那里的算法。从页面:

程序员:本文档中的 JavaScript 源代码可以无限制地复制和重复使用。

链接好像坏了!
2021-04-21 15:47:20
这是我见过的最干净的 Javascript 版本。
2021-04-22 15:47:20
正如@monkut 的回答中提到的,python.org 包括来自 python 的 GDAL 的 Swig 绑定:pypi.python.org/pypi/GDAL
2021-05-15 15:47:20

根据此页面,proj4js 支持 UTM。

http://trac.osgeo.org/proj4js/wiki/UserGuide#Supportedprojectionclasses

您可能还想看看GDALgdal 库具有出色的 python 支持,但如果您只进行投影转换,它可能有点过分。

+1 PROJ4 几乎可以做您梦寐以求的任何事情,所以如果 proj4js 是一个真正的端口,它也可以做到。
2021-04-22 15:47:20
+1 proj4s 是要走的路。没有必要重新发明轮子.. proj4s 从可以添加任何投影的配置文件中读取 - 将引用输入到spatialreference.org 中以获取 proj4js 字符串,例如spatialreference.org/ref/epsg/4326/proj4js
2021-05-07 15:47:20

我也是新手,最近一直在研究这个主题。

这是我使用 python gdal pacakge找到的方法osr包包含在 gdal 中)。gdal 包非常强大,但文档可能会更好。

这源自此处的讨论:http : //www.mail-archive.com/gdal-dev@lists.osgeo.org/msg12398.html

import osr

def transform_utm_to_wgs84(easting, northing, zone):
    utm_coordinate_system = osr.SpatialReference()
    utm_coordinate_system.SetWellKnownGeogCS("WGS84") # Set geographic coordinate system to handle lat/lon
    is_northern = northing > 0    
    utm_coordinate_system.SetUTM(zone, is_northern)

    wgs84_coordinate_system = utm_coordinate_system.CloneGeogCS() # Clone ONLY the geographic coordinate system 

    # create transform component
    utm_to_wgs84_transform = osr.CoordinateTransformation(utm_coordinate_system, wgs84_coordinate_system) # (<from>, <to>)
    return utm_to_wgs84_transform.TransformPoint(easting, northing, 0) # returns lon, lat, altitude

这是从 wgs84(大多数 gps 单位报告的内容)中的经纬度转换为 utm 的方法:

def transform_wgs84_to_utm(lon, lat):    
    def get_utm_zone(longitude):
        return (int(1+(longitude+180.0)/6.0))

    def is_northern(latitude):
        """
        Determines if given latitude is a northern for UTM
        """
        if (latitude < 0.0):
            return 0
        else:
            return 1

    utm_coordinate_system = osr.SpatialReference()
    utm_coordinate_system.SetWellKnownGeogCS("WGS84") # Set geographic coordinate system to handle lat/lon  
    utm_coordinate_system.SetUTM(get_utm_zone(lon), is_northern(lat))

    wgs84_coordinate_system = utm_coordinate_system.CloneGeogCS() # Clone ONLY the geographic coordinate system 

    # create transform component
    wgs84_to_utm_transform = osr.CoordinateTransformation(wgs84_coordinate_system, utm_coordinate_system) # (<from>, <to>)
    return wgs84_to_utm_transform.TransformPoint(lon, lat, 0) # returns easting, northing, altitude    

我还发现,如果您已经安装了 django/gdal 并且您知道您正在处理的 UTM 区域EPSG代码,您可以只使用Point() transform()方法。

from django.contrib.gis.geos import Point
utm2epsg = {"54N": 3185, ...}
p = Point(lon, lat, srid=4326) # 4326 = WGS84 epsg code
p.transform(utm2epsg["54N"])

Staale 答案的 Javascript 版本

function utmToLatLng(zone, easting, northing, northernHemisphere){
        if (!northernHemisphere){
            northing = 10000000 - northing;
        }

        var a = 6378137;
        var e = 0.081819191;
        var e1sq = 0.006739497;
        var k0 = 0.9996;

        var arc = northing / k0;
        var mu = arc / (a * (1 - Math.pow(e, 2) / 4.0 - 3 * Math.pow(e, 4) / 64.0 - 5 * Math.pow(e, 6) / 256.0));

        var ei = (1 - Math.pow((1 - e * e), (1 / 2.0))) / (1 + Math.pow((1 - e * e), (1 / 2.0)));

        var ca = 3 * ei / 2 - 27 * Math.pow(ei, 3) / 32.0;

        var cb = 21 * Math.pow(ei, 2) / 16 - 55 * Math.pow(ei, 4) / 32;
        var cc = 151 * Math.pow(ei, 3) / 96;
        var cd = 1097 * Math.pow(ei, 4) / 512;
        var phi1 = mu + ca * Math.sin(2 * mu) + cb * Math.sin(4 * mu) + cc * Math.sin(6 * mu) + cd * Math.sin(8 * mu);

        var n0 = a / Math.pow((1 - Math.pow((e * Math.sin(phi1)), 2)), (1 / 2.0));

        var r0 = a * (1 - e * e) / Math.pow((1 - Math.pow((e * Math.sin(phi1)), 2)), (3 / 2.0));
        var fact1 = n0 * Math.tan(phi1) / r0;

        var _a1 = 500000 - easting;
        var dd0 = _a1 / (n0 * k0);
        var fact2 = dd0 * dd0 / 2;

        var t0 = Math.pow(Math.tan(phi1), 2);
        var Q0 = e1sq * Math.pow(Math.cos(phi1), 2);
        var fact3 = (5 + 3 * t0 + 10 * Q0 - 4 * Q0 * Q0 - 9 * e1sq) * Math.pow(dd0, 4) / 24;

        var fact4 = (61 + 90 * t0 + 298 * Q0 + 45 * t0 * t0 - 252 * e1sq - 3 * Q0 * Q0) * Math.pow(dd0, 6) / 720;

        var lof1 = _a1 / (n0 * k0);
        var lof2 = (1 + 2 * t0 + Q0) * Math.pow(dd0, 3) / 6.0;
        var lof3 = (5 - 2 * Q0 + 28 * t0 - 3 * Math.pow(Q0, 2) + 8 * e1sq + 24 * Math.pow(t0, 2)) * Math.pow(dd0, 5) / 120;
        var _a2 = (lof1 - lof2 + lof3) / Math.cos(phi1);
        var _a3 = _a2 * 180 / Math.PI;

        var latitude = 180 * (phi1 - fact1 * (fact2 + fact3 + fact4)) / Math.PI;

        if (!northernHemisphere){
          latitude = -latitude;
        }

        var longitude = ((zone > 0) && (6 * zone - 183.0) || 3.0) - _a3;

        var obj = {
              latitude : latitude,
              longitude: longitude
        };


        return obj;
      }