解锁空间距离计算的多种方式-含前端、空间数据库、后端

前言

一、空间数据库求解

1、PostGIS实现

二、GIS前端组件求解

1、Leaflet.js距离测算

2、Turf.js前端计算

三、后台距离计算生成

1、欧式距离

2、Haversice球面距离

3、GeoTools距离计算

4、Gdal距离生成

5、geodesy距离计算

四、成果与生成对比

1、Java不同生成方法对比

2、各种生成方式对比

五、总结

前言

空间距离是指物体在空间中的位置之间的距离，通常用来描述物体之间的相对位置关系。在日常生活中，我们经常使用距离来描述物体的位置关系，比如在行驶中使用路程来描述两个地点之间的距离，或者在导航中使用地图上的距离来指引行驶方向。在物理学和数学中，距离是一个重要的概念，它被用来描述空间中的位置关系，衡量物体之间的远近。空间距离的研究对于理解物体的位置关系、运动轨迹、引力场等具有重要的意义。下面是在某地图软件中，使用测距的方式直接量算的长度示例。

在地理世界中，在导航系统中，物流系统中，旅游景点等应用。距离这个概念更是非常常见的，比如在导航系统中，从家到某商场的距离；在快递物流行业中，我们从浙江义乌的小商品城购买一件物品，然后邮寄到家；还可能我们计划出门旅游，那么我们想知道家到目的地的实际距离有多远等等场景，都是非常具有代表性的距离求解和应用场景。

伴随着生活的不断进步，在我们的生活中，我们会应用更多距离计算。作为一位GIS开发者，距离是空间分析中最简单的一种类型。本文将详细的讲解各种不同位置距离计算方法，首先讲解使用空间数据库的直接求解办法，其次介绍在前端组件如Leaflet.js、Turf.js等组件中进行求解的办法，然后基于Java语言讲解后端的详细计算方法，包括Java直接求解、GeoTools计算、geodesy距离计算、欧式距离计算、GDAL空间计算等不同的方法，最后对比不同的计算方法得到的结果，为大家对距离的计算有更多的掌握和了解。

一、空间数据库求解

在进行空间距离求解时，我们这里假定一种场景。比如以共享单车为例，需要计算从起始点到结束点的距离，这里涉及到两个信息，也就是已知两个经纬度坐标点，即： $p_{1}\left (x_{1},y_{1} \right )$ 和 $p_{2}\left (x_{2},y_{2} \right )$ ，在给定两个包含经纬度坐标的点后，我们来求解其距离。我们计算从某公园的一角到地铁站的距离，首先我们在地图上将这个点标记出来，得到这两点的实际经纬度。实现方法是在Leaflet中绑定点击事件，然后在点击事件中输出具体的坐标。

在地图上标记出两个位置，得到需要计算的两个点的具体经纬度信息： $p_{1}\left (112.867699, 28.194674 \right )$ 和 $p_{2}\left (112.876024, 28.189038 \right )$ 。 $p_{1}$ 表示起点， $p_{2}$ 表示目的地。后续的计算均以这两个坐标点为起始位置。

1、PostGIS实现

为了简单起见，这里的空间数据库实现，我们采用PostGIS来进行开发数据库。当然在实际开发过程中，除了可以选择PostGIS数据库，还可以选择其它的空间数据库，比如MySQL或者Oracle SDE等。在PostGIS数据库中，我们可以使用ST_Distance()和ST_GeomFromText()进行查询，首先使用ST_GeomFromText()将字符串转换为Geography，然后再使用ST_Distance()来进行记录求解，请注意，在进行求解时，务必保证安装PostGIS扩展。

SELECT
	ST_Distance ( 
    ST_GeomFromText ( 'POINT(112.867699 28.194674)', 4326 ) :: geography, 
	ST_GeomFromText ( 'POINT(112.876024 28.189038)', 4326 ) :: geography 
) AS distance;

这里选择将geometry转换为geography字段是为了将输出的结果转换为常见的距离单位，米。在postgis中，会自动根据不同的投影类型来得到距离结果，如果用4326的参考系，得到的结果将是度。在PostGIS中使用空间数据库的方式的消耗时间大约是：3毫秒，两点间的测算距离是1028.7251米。速度也是非常快的。

二、GIS前端组件求解

在介绍完空间数据库的求解模式后。接下来我们介绍一下基于WebGIS的前端组件求解方式。这种解决方式比较适用于以前端为主的应用中，不需要数据库的支持。常见的WebGIS前端组件也有许多，二维常见的有Leaflet、OpenLayers等，三维有Cesium等。这里我们已Leaflet为例，同时讲解在Turf.js中是如何实现的。

1、Leaflet.js距离测算

首先来看一下再Leaflet.js中进行距离的计算，得到的结果是多少。首先我们需要在页面中引用到Leaflet.js，可以是本地的离线资源，也可以是在线资源。可以直接使用Leaflet.js的相关API来进行。关键代码如下：

let latlng1 = L.latLng(28.194674, 112.867699);
// 第二个点
let latlng2 = L.latLng(28.189038, 112.876024);
// 计算两点之间的距离
let distance = latlng1.distanceTo(latlng2);

使用Leaflet的距离求解方法得到的距离是：1028.79098米，耗时5毫秒。

2、Turf.js前端计算

Turf.js也是一款在WebGIS系统开发当中很常见的前端计算框架，在之前的一些博客中曾经做了比较详细的介绍，关于turf.js，它就像一个宝藏一样，值得我们去学习研究，发现更多的实用功能。turf可以和Leaflet或者openLayers，cesium等配合使用。这里重点讲解如何使用Turf.js来进行距离运算。关键代码如下：

var from = turf.point([112.867699,28.194674]);
var to = turf.point([112.876024,28.189038]);
var options = { units: "meters" };//设置单位为米
var distance = turf.distance(from, to, options);//求解计算

使用Turf.js的距离求解方法得到的距离是：1028.7924米，耗时4毫秒。

三、后台距离计算生成

在了解空间数据库生成和WebGIS前端组件生成方式之后，我们来看一下如何使用后台服务的生成方式。这里以Java开发语言为例，主要介绍如何使用Java来进行不同的距离进行生成。包括但不限于欧式距离、球面距离、GeoTools距离生成、Gdal距离生成、Geodesy距离生成等。

1、欧式距离

欧式距离（Euclidean distance）是最直接距离度量方式，用于计算两个点在欧几里得空间中的直线距离，如果有两个点： $p_{1}\left ( x_{1},y_{1} \right )$ 和 $p_{2}\left ( x_{2},y_{2} \right )$ 在二维空间中，它们之间的欧式距离d 可以通过下面的公式来计算：

$d \left (p_{1},p_{2} \right )=\sqrt{\left ( x_{2}-x_{1} \right )^{2}+\left ( y_{2}-y_{1} \right )^{2}}$

根据上面的公式定义，可以直接使用Java原始的计算方式来计算，关键代码如下所示：

/**
* 计算两点的欧式距离
* @param p0  起始点
* @param p1 结束点
* @return
*/
public static double calculateEuclideanDistance(Coordinate p0, Coordinate p1) {
    Long start = System.currentTimeMillis();
    double result = Math.sqrt(Math.pow(p1.x - p0.x, 2) + Math.pow(p1.y - p0.y, 2));
    Long end= System.currentTimeMillis();
    System.out.println("计算耗时：" + ( end - start)  + "毫秒");
    return result;
 }

2、Haversice球面距离

Haversice公式，通过两点计算大圆距离，（即球面上两点的最短距离），Haversice的数学表达：设两点的经纬度分别为（lat1，lon1）和（lat2，lon2），其中经度和纬度以弧度为单位。地球的平均半径R约为6371公里。Haversice的计算步骤为：

首先将经纬度从度数转换为弧度（考虑曲率）。lat1,lon1,lat2,lon2 ->rad。然后计算经度和纬度的差值， $\Delta lat=lat_{2}-lat_{1}$ 以及 $\Delta lon=lon_{2}-lon_{1}$ ,最后使用下面的公式计算距离：

$a=\sin ^{2}\left ( \frac{\Delta lat}{2} \right )+\cos \left ( lat_{1} \right )+\cos \left ( lat_{2} \right )+\sin ^{2}\left ( \frac{\Delta lon}{2} \right )$

$c=2\times \arctan 2\left ( \sqrt{a},\sqrt{1-a} \right )$

d=R * c

在上面的公式中，d表示两点间的距离。将上述公式转换为代码后如下所示：

/**
 * 计算两个位置的球面距离
 * @param p0 起始点
 * @param p1 目标点
 * @return
 */
public static double calculateHaversineDistance(Coordinate p0, Coordinate p1) {
   Long start = System.currentTimeMillis();
   double radLat1 = Math.toRadians(p0.x);
   double radLon1 = Math.toRadians(p0.y);
   double radLat2 = Math.toRadians(p1.x);
   double radLon2 = Math.toRadians(p1.y);

   double deltaLat = radLat2 - radLat1;
   double deltaLon = radLon2 - radLon1;

   double a = Math.pow(Math.sin(deltaLat / 2), 2) +  Math.cos(radLat1) * Math.cos(radLat2) *
                   Math.pow(Math.sin(deltaLon / 2), 2);
   double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
   double distance = EARTH_RADIUS * c;
   Long end= System.currentTimeMillis();
   System.out.println("计算耗时：" + ( end - start)  + "毫秒");
   return distance;
}

3、GeoTools距离计算

我们可以直接利用GeoTools来直接进行距离的计算，它以API的方式直接提供调用，这里我们不深入核心的计算方法，具体大家可以深入源码去看。示例的关键代码如下：

/**
* 使用Geotools来求解距离使用WGS84参考系
* @param p0 起始点
* @param p1 目标点
* @return
*/
public static double geoToolsDistanceWithWgs84Crs(Coordinate p0, Coordinate p1) {
	Long start = System.currentTimeMillis();
	// 84坐标系构造GeodeticCalculator
    GeodeticCalculator geodeticCalculator = new GeodeticCalculator(DefaultGeographicCRS.WGS84);
    // 起点经纬度
    geodeticCalculator.setStartingGeographicPoint(p0.x, p0.y);
    // 末点经纬度
    geodeticCalculator.setDestinationGeographicPoint(p1.x, p1.y);
    // 计算距离，单位：米
    double result = geodeticCalculator.getOrthodromicDistance();
    Long end= System.currentTimeMillis();
    System.out.println("计算耗时：" + ( end - start) + "毫秒");
    return result;
 }

4、Gdal距离生成

与GeoTools类似，在Java开发当中，除了可以使用GeoTools中进行距离计算。Gdal也是同样可以进行距离计算的。

 /**
  * 使用gdal来计算两个点距离
  * @param p0 起始点
  * @param p1 目标点
  * @return
  */
 public static double calculateGdalDistance(Coordinate p0, Coordinate p1) {
    Long start = System.currentTimeMillis();
    // 注册所有的GDAL驱动
    gdal.AllRegister();
    // 创建两个点的几何对象
    Geometry point1 = ogr.CreateGeometryFromWkt("POINT ("+ p0.y+" "+p0.x+")");
    Geometry point2 = ogr.CreateGeometryFromWkt("POINT ("+ p1.y+" "+p1.x+")");
        
    // 创建空间参考对象
    SpatialReference srs = new SpatialReference();
    // 导入WGS84坐标系的EPSG代码
    srs.ImportFromEPSG(4326);
    // 为几何对象设置空间参考
    point1.AssignSpatialReference(srs);
    point2.AssignSpatialReference(srs);

    // 计算两个点之间的距离
    double distance = point1.Distance(point2);
    // 输出距离，注意GDAL返回的距离单位可能与坐标系单位不同，需要根据实际情况转换
    System.out.println("Distance: " + distance);
    Long end= System.currentTimeMillis();
    System.out.println("计算耗时：" + ( end - start)  + "毫秒");
    return distance;
}

5、geodesy距离计算

与Gdal与GeoTools不同，geodesy是一个轻量级的地理分析和处理框架。因此我们可以可以使用geodesy来进行距离的计算，关键计算过程如下所示：

/**
 * 使用geodesy来计算两个点距离
 * @param p0 起始点
 * @param p1 目标点
 * @return
 */
 public static double calculateGeodesyDistance(Coordinate p0, Coordinate p1) {
    Long start = System.currentTimeMillis();
    //创建GeodeticCalculator，调用计算方法，传入坐标系、经纬度用于计算距离
    GeodeticCurve geoCurve = new org.gavaghan.geodesy.GeodeticCalculator().calculateGeodeticCurve(
                Ellipsoid.WGS84,
                new GlobalCoordinates(p0.x, p0.y),
                new GlobalCoordinates(p1.x, p1.y)
    );
    double result = geoCurve.getEllipsoidalDistance();
    Long end= System.currentTimeMillis();
    System.out.println("计算耗时：" + ( end - start)  + "毫秒");
    return result;
}

以上就是5种使用java生成距离的方式，通过以上的代码均可生成距离。

四、成果与生成对比

这里将对生成的距离成果与不同的生成方式进行对比，方便大家对结果和生成过程有一个具体的认识。

1、Java不同生成方法对比

我们以Junit代码为例，分别调用上述的代码，用来计算两个点的距离和生成时间。代码如下：

@Test
public void testGeoToolsDis() {
	Coordinate coord1 = new Coordinate(112.867699, 28.194674); //出发的位置
    Coordinate coord2 = new Coordinate(112.876024, 28.189038); // 到达的经纬度 
    double distance = geoToolsDistanceWithDefaultCrs(coord1, coord2);
    System.out.println("geotools Distance no crs: " + distance + " meters");
    System.out.println("********************************************");
    distance = geoToolsDistanceWithWgs84Crs(coord1,coord2);
    System.out.println("geotools Distance crs: " + distance + " meters");
    System.out.println("********************************************");
    double distance1 = calculateEuclideanDistance(coord1,coord2);
    System.out.println("欧式距离(度)：" + distance1);
    distance = distance1 * EARTH_RADIUS * Math.PI / 180 * 1000;
    System.out.println("欧式距离(米）：" + distance + " meters");
    System.out.println("********************************************");
    distance = calculateHaversineDistance(coord1,coord2) * 1000;
    System.out.println("球面距离：Haversine Distance: " + distance + " meters");
    System.out.println("********************************************");
    distance = calculateGeodesyDistance(coord1,coord2);
    System.out.println("Geodesy Distance: " + distance + " meters");
    System.out.println("********************************************");
    distance = calculateGdalDistance(coord1,coord2);
    distance = distance * EARTH_RADIUS * Math.PI / 180 * 1000;
    System.out.println("Gdal Distance: " + distance + " meters");
}

可以控制台看到以下的输出：

计算耗时：230毫秒
geotools Distance no crs: 1028.725087712003 meters
********************************************
计算耗时：1790毫秒
geotools Distance crs: 1028.725087712003 meters
********************************************
计算耗时：0毫秒
欧式距离(度)：0.010053363665956643
欧式距离(米）：1119.1353240538915 meters
********************************************
计算耗时：0毫秒
球面距离：Haversine Distance: 958.2800850973537 meters
********************************************
计算耗时：5毫秒
Geodesy Distance: 960.1006447758463 meters
********************************************
Distance: 0.010053363665956643
计算耗时：336毫秒
Gdal Distance: 1119.1353240538915 meters

请注意，在上面的计算过程中，使用欧式距离是，首先生成出来的是度，而不是我们的预期米。因此需要进行转换。转换代码可以参考上面的示例代码。

2、各种生成方式对比

这里将提供两种数据结果对比，第一种是提供表格的形式，第二种是提供图表的形式。如下表所示：

序号	生成方式	生成距离（米）	耗时（毫秒）
1	PostGIS数据库	1028.72508	3
2	Leaflet.js	1028.79098	5
3	Turf.js	1028.7924	4
4	GeoTools no Crs	1028.725087	230
5	GeoTools with Crs	1028.725087	1790
6	欧式距离	1119.13532	接近0
7	Haversine距离	958.28008	接近0
8	Geodesy	960.10064	5
9	Gdal 距离	1119.13532	336

将上面的各种距离生成结果使用图表的形式来展示，如下所示：

五、总结

以上就是本文的主要内容，本文将详细的讲解各种不同位置距离计算方法，首先讲解使用空间数据库的直接求解办法，其次介绍在前端组件如Leaflet.js、Turf.js等组件中进行求解的办法，然后基于Java语言讲解后端的详细计算方法，包括Java直接求解、GeoTools计算、geodesy距离计算、欧式距离计算、GDAL空间计算等不同的方法，最后对比不同的计算方法得到的结果，为大家对距离的计算有更多的掌握和了解。行文仓促，定有许多不足之处，如有不足，还请各位专家博主在评论区留言指出，不胜感激。