屏幕xy相对坐标转经纬度坐标

2,779 阅读7分钟

今天读了一下《地理坐标(经纬度坐标)和屏幕坐标(xy坐标)间的转换》,觉得内容很好,特摘录和附上自己的解释。

背景

在我们的屏幕上,有一张地图,这张地图经过缩放、平移、旋转,最终地理坐标和屏幕坐标的关系大致如下图所示:

这种关系要怎么描述呢?我们可以假设地图是一张纸,而屏幕是一堵墙。只要我们有两个图钉,我们就能把纸定在墙上。我们把这两个点称为锚点。锚点在屏幕坐标系上的坐标是(x1,y1)(x2,y2),对应在地理坐标系上的坐标是(lon1,lat1)(lon2,lat2)

那现在的问题就变成了,已知两个锚点的坐标,

(1)地理坐标转屏幕坐标:已知任意一点的地理坐标(lon,lat),求它在屏幕上的坐标(x,y)

(2)屏幕坐标转地理坐标:已知任意一点的屏幕坐标(x,y),求它的经纬度坐标(lon,lat)

lon 表示 longitude (纬度)
lat 表示 latitude (经度)
x,y 表示屏幕上的只是一个任意位置的坐标点

转换算法

1、铺垫内容-地理坐标平面化

在一个小范围内(例如是方圆几公里内),我们可以假设地面是平的,而不是弯的。如果经纬度都用弧度表示,那么1纬度对应的长度是:

1纬度长度=Rlat1 纬度长度 = R \cdot lat

注意:此处的 lat 指的不是我们常用的纬度度数,而是纬度弧度。半径乘以弧度才能得到弧长! 度数360=弧度2π\frac{度数}{360}=\frac{弧度}{2π} 故 弧度 = 度数 π180\cdot \frac{ \pi }{180}

其中 RR 是地球半径。 而相同经度间的距离会随着纬度的增加而减少,

image.png

latlat 这一纬度下,1经度对应的长度是:

1经度长度=Rloncos(lat)1 经度长度= R \cdot lon \cdot \cos{(lat )}

上式中的 Rcos(lat)R \cdot \cos{(lat )} 就是上图中的某一纬度下的切面的半径,然后在半径乘以弧度。

那么,(lon,lat)这个坐标平面化后的坐标就是: (Rloncos(lat),Rlat)(R \cdot lon \cdot \cos{(lat)},R \cdot lat)

由于我们是假定在小范围内纬度是不变的,故此时若我们有两个经度不同的点,求经度长度的差时,应该遵循:

dlenlon=Rlon2cos(lat)Rlon1cos(lat)=R(lon2lon1)cos(lat)dlen_{lon} = R \cdot lon_2 \cdot cos (lat) - R \cdot lon_1 \cdot cos (lat) = R(lon_2 - lon_1 )\cdot cos (lat)

其实上式中的求纬度长度差时用到的latlat,使用lat1lat_1或者是lat2lat_2都对,因为我们在计算经度长度时视作一致。

如果求一个纬度长度差与经度长度差比值的话:

dlenlatdlenlon=Rlat2Rlat1R(lon2lon1)cos(lat)\frac{dlen_{lat}}{dlen_{lon}} = \frac{R \cdot lat_2-R \cdot lat_1}{R(lon_2 - lon_1 )\cdot cos (lat)}

好,我们上下同时除以RR,可得

=lat2lat1(lon2lon1)cos(lat)= \frac{lat_2- lat_1}{(lon_2 - lon_1 )\cdot cos (lat)}

又因为前面提到过,我们计算经度长度时视纬度为不变的,故我们随意地,把 lat2lat_2 代给 latlat .

=lat2lat1(lon2lon1)cos(lat2)= \frac{lat_2- lat_1}{(lon_2 - lon_1 )\cdot cos (lat_2)}

这里需要注意的是:lat2lat1(lon2lon1)\frac{lat_2- lat_1}{(lon_2 - lon_1 )} 是个比值,故用度数比和弧度比都一致。但是,cos(lat2)cos (lat_2) 处只能是弧度,原文的代码中有特别指出。

此二者的比值可以方便我们:

在得到经度时根据比值关系获得纬度

在得到纬度时根据比值关系获得经度

换到解方程组里就是y通过某种关系变成了能用x表示的式子,方程中只剩下一个未知数后,便可以求解。此刻,我们已经得到经纬度的代数关系,下面就是去根据向量法来把问题拉回到一个点上的经纬度以及xy之间的关系。

2、核心内容-向量法

由已知点和未知点组成两组向量:

由于坐标系转换是线性变换,所以两组向量有以下特性:

  • (1)两向量在不同的坐标系中的长度比是相同的。

  • (2)两向量在不同的坐标系中的夹角是相同的。

根据上面两个特性,我们可列出方程组:

下面提到的 d ,应该只是简单代指 difference (差值)

设向量1为(dx1,dy1)(dx_1,dy_1)(dlon1,dlat1)(dlon_1,dlat_1),向量2为 (dx2,dy2)(dlon2,dlat2)(dx_2,dy_2),(dlon_2,dlat_2)

其中

dx1=x2x1dx_1=x_2-x_1dy1=y2y1dy_1=y_2-y_1dlon1=lon2lon1dlon_1=lon_2-lon_1dlat1=lat2lat1dlat_1=lat_2-lat_1

dx2=xx1dx_2=x-x_1dy2=yy1dy_2=y-y_1dlon2=lonlon1dlon_2=lon-lon_1dlat2=latlat1dlat_2=lat-lat_1

然后

k1=norm(dx1,dy1)k_1=norm(dx_1,dy_1)

k2=norm(dlon1,dlat1)k_2=norm(dlon_1,dlat_1)

k3=norm(dx2,dy2)k_3=norm(dx_2,dy_2)

k4=norm(dlon2,dlat2)k_4=norm(dlon_2,dlat_2)

这里的 norm 函数我查了一下,就是求范数,根据《np.linalg.norm(求范数)》 一文的介绍,向量的范数有计算公式,默认是求二范数,二范数就是向量的模长

下面我们给出方程组:

  1. 根据前文中提到的 两向量在不同的坐标系中的长度比是相同的 可得到
k1k2=k3k4\frac{k_1}{k_2} = \frac{k_3}{k_4}

本式当然可以写成:

dx12+dy12dlon12+dlat12=dx22+dy22dlon22+dlat22\frac{\sqrt{dx_1^2 + dy_1^2}}{\sqrt{dlon_1^2+dlat_1^2}} = \frac{\sqrt{dx_2^2 + dy_2^2}}{\sqrt{dlon_2^2+dlat_2^2}}
  1. 根据前文中提到的 两向量在不同的坐标系中的夹角是相同的

我们可得出以下四个式子:

dx1k1=dx2k3(通过x求夹角cos值)① \frac{dx_1}{k_1} = \frac{dx_2}{k_3} (通过x求夹角 cos 值)
dy1k1=dy2k3(通过y求夹角sin值)② \frac{dy_1}{k_1} = \frac{dy_2}{k_3} (通过y求夹角 sin 值)
dlon1k2=dlon2k4(通过经度求夹角cos值)③ \frac{dlon_1}{k_2} = \frac{dlon_2}{k_4} (通过经度求夹角 cos 值)
dlat1k2=dlat2k4(通过纬度求夹角sin值)④ \frac{dlat_1}{k_2} = \frac{dlat_2}{k_4} (通过纬度求夹角 sin 值)

由 ① 乘以 ③ 可得:

dx1dlon1k1k2=dx2dlon2k3k4⑤ \frac{dx_1 \cdot dlon_1 }{k1 \cdot k_2} = \frac{dx_2 \cdot dlon_2}{k_3 \cdot k_4}

由 ② 乘以 ④ 可得:

dy1dlat1k1k2=dy2dlat2k3k4⑥ \frac{dy_1 \cdot dlat_1 }{k1 \cdot k_2} = \frac{dy_2 \cdot dlat_2}{k_3 \cdot k_4}

将 ⑤ 与 ⑥ 相加,即可得原文中提到的复合在一起的式子:

dx1dlon1+dy1dlat1k1k2=dx2dlon2+dy2dlat2k3k4\frac{dx_1 \cdot dlon_1 + dy_1 \cdot dlat_1 }{k1 \cdot k_2} = \frac{dx_2 \cdot dlon_2 + dy_2 \cdot dlat_2}{k_3 \cdot k_4}

故最后的方程组为:

{k1k2=k3k4dx1dlon1+dy1dlat1k1k2=dx2dlon2+dy2dlat2k3k4\begin{array}{l} \left\{\begin{matrix} \frac{k_1}{k_2} = \frac{k_3}{k_4} \\ \frac{dx_1 \cdot dlon_1 + dy_1 \cdot dlat_1 }{k1 \cdot k_2} = \frac{dx_2 \cdot dlon_2 + dy_2 \cdot dlat_2}{k_3 \cdot k_4} \end{matrix}\right. \end{array}

我再贴张用LaTeX公式编辑器绘制出的大图:

请添加图片描述

通过解上面的方程组,我们就能得到未知和屏幕坐标或未知的地理坐标。

为什么呢?别看推导过程中代数符号(变量名)很多,其实不确定的值只有:

(x,y)(lon,lat)(x,y)和(lon,lat)

故,接下来就简单了:

  • (x,y)(x,y)带入方程组就可以得到(lon,lat)(lon,lat)

  • 反过来,将(lon,lat)(lon,lat)带入方程组就可以得到(x,y)(x,y)


这里需要注意的是:

处于分母位置的 k1k_1k2k_2k3k_3k4k_4 都不能为零

转换为编程中的除零异常就是需要注意的是:

  • 参考的两点不能完全相等
  • 进行计算的(x,y)(x,y)(lon,lat)(lon,lat)不能与第一个参考点相同

对于第二句话,由于类库使用者调用本函数时,不会注意传参顺序,故你可要求他要计算的点不应该为两个参考点中任意一个。 实际工程中,对于这个问题,我们判断一下 dx1,dy1dx_1,dy_1 ,即当前点是否和参考点重合,如果重合,直接返回参考点。

屏幕坐标转地理坐标(Java实现)

由于原文中已给出了详尽的C#代码,这里我只给出自己的关于屏幕坐标转地理坐标的Java实现。由于代码实现解方程组我还不熟悉,只能按照原有文章中的代码逻辑直接给出到大家,待后续有能力了,我会把变量都和代数中的符号保持一致。

1.首先定义包含(x,y)(x,y)(lon,lat)(lon,lat) 的数据存储的类

public class XYLonLatModle {

    private double x;

    private double y;

    private double longitude;

    private double latitude;

    public double getX() {
        return x;
    }

    public void setX(double x) {
        this.x = x;
    }

    public double getY() {
        return y;
    }

    public void setY(double y) {
        this.y = y;
    }

    public double getLongitude() {
        return longitude;
    }

    public void setLongitude(double longitude) {
        this.longitude = longitude;
    }

    public double getLatitude() {
        return latitude;
    }

    public void setLatitude(double latitude) {
        this.latitude = latitude;
    }

    public String getWgsString() {
        return latitude + "," + longitude;
    }
}

然后为了构建这个 XYLonLatModle 更加方便,我们再写一个 Builder 类。

public final class XYLonLatModleBuilder {
    private double x;
    private double y;
    private double longitude;
    private double latitude;

    private XYLonLatModleBuilder() {
    }

    public static XYLonLatModleBuilder aXYLonLatModle() {
        return new XYLonLatModleBuilder();
    }

    public XYLonLatModleBuilder withX(double x) {
        this.x = x;
        return this;
    }

    public XYLonLatModleBuilder withY(double y) {
        this.y = y;
        return this;
    }

    public XYLonLatModleBuilder withOriginalXY() {
        this.x = 0;
        this.y = 0;
        return this;
    }

    public XYLonLatModleBuilder withXyString(String xyString){
        String[] strings = xyString.split(",");
        this.x = Double.parseDouble(strings[0]);
        this.y = Double.parseDouble(strings[1]);
        return this;
    }

    public XYLonLatModleBuilder withLongitude(double longitude) {
        this.longitude = longitude;
        return this;
    }

    public XYLonLatModleBuilder withLatitude(double latitude) {
        this.latitude = latitude;
        return this;
    }

    public XYLonLatModleBuilder withLatitudeLongitudeString(String latLngString) {
        String[] strings = latLngString.split(",");
        this.latitude = Double.parseDouble(strings[0]);
        this.longitude = Double.parseDouble(strings[1]);
        return this;
    }

    public XYLonLatModle build() {
         XYLonLatModle xYLonLatModle = new  XYLonLatModle();
        xYLonLatModle.setX(x);
        xYLonLatModle.setY(y);
        xYLonLatModle.setLongitude(longitude);
        xYLonLatModle.setLatitude(latitude);
        return xYLonLatModle;
    }
}

2.创建工具类根据两个参照点和当前的(x,y)(x,y) 计算 (lon,lat)(lon,lat)

public class XYLongitudeLatitudeConvertUtil {

    public static double norm(double a, double b) {
        return Math.sqrt(a * a + b * b);
    }

    /**
     * x y 坐标转 经纬度
     *
     * @param x
     * @param y
     * @param reference1
     * @param reference2
     * @return
     */
    public static XYLonLatModle vectorMapXyToLatLng(double x, double y, XYLonLatModle reference1, XYLonLatModle reference2) {

        double x1 = reference1.getX();
        double y1 = reference1.getY();
        double lon1 = reference1.getLongitude();
        double lat1 = reference1.getLatitude();

        double x2 = reference2.getX();
        double y2 = reference2.getY();
        double lon2 = reference2.getLongitude();
        double lat2 = reference2.getLatitude();
        

        if (x - x1 == 0 && y - y1 == 0) {
            // 实时点 与第一参考点重合,直接返回参考点
            return reference1;
        }
        if (x - x2 == 0 && y - y2 == 0) {
            // 实时点 与第二参考点重合,直接返回参考点
            return reference2;
        }


        double lon_cos = Math.cos(lat2 * Math.PI / 180);

        double m = (lon2 - lon1) * lon_cos;
        double n = (lat2 - lat1);

        double dx1 = x2 - x1;
        double dy1 = y2 - y1;
        double dx2 = x - x1;
        double dy2 = y - y1;

        double a = (dx2 * dx2 + dy2 * dy2) * (m * m + n * n) / (dx1 * dx1 + dy1 * dy1);

        double norm = norm(m, n);
        double v = dx1 * dx2 + dy2 * dy1;
        double sqrt = Math.sqrt(a);
        double b = v * norm * sqrt / (norm(dx1, dy1) * norm(dx2, dy2));

        double c = Math.sqrt(b * b * n * n - (m * m + n * n) * (b * b - a * m * m));

        double q1 = (b * n + c) / (m * m + n * n);
        double q2 = (b * n - c) / (m * m + n * n);

        double p1 = (b - q1 * n) / m;
        double p2 = (b - q2 * n) / m;

        double lon_1 = p1 / lon_cos + lon1;
        double lat_1 = q1 + lat1;
        double lon_2 = p2 / lon_cos + lon1;
        double lat_2 = q2 + lat1;

        double judge1 = (lon_1 - lon1) * (lat2 - lat1) - (lat_1 - lat1) * (lon2 - lon1);
        double judge2 = (lon_2 - lon1) * (lat2 - lat1) - (lat_2 - lat1) * (lon2 - lon1);
        double judge = (x - x1) * (y2 - y1) - (y - y1) * (x2 - x1);

        double lon = 0;
        double lat = 0;
        if (judge * judge1 < 0) {
            lon = lon_1;
            lat = lat_1;
        } else {
            lon = lon_2;
            lat = lat_2;
        }

        XYLonLatModle xyLonLatModle = XYLonLatModleBuilder.aXYLonLatModle().withX(x).withY(y).withLongitude(lon).withLatitude(lat).build();
        return xyLonLatModle;
    }

}

本方法未用到任何jar包依赖,直接使用的Java基本语法,故不需要给出import

最简易的调用就是用屏幕中的零点和一个较大的点来作参考:

XYLonLatModle originalPoint = XYLonLatModleBuilder.aXYLonLatModle()
        .withOriginalXY()
        .withLatitudeLongitudeString("20.156,52.369")
        .build();
XYLonLatModle anotherPoint = XYLonLatModleBuilder.aXYLonLatModle()
        .withXyString("5,5")
        .withLatitudeLongitudeString("20.157,52.370")
        .build();
XYLonLatModle xyLonLatModle = XYLongitudeLatitudeConvertUtil.vectorMapXyToLatLng(x, y, originalPoint, anotherPoint);