计算船只到岸或海岸线的距离

Question

对于一个200M GPS（lon，lat）船只坐标的数据集，我想计算到最近的陆地或海岸线的近似距离，作为名为distance_to_shore的函数，该函数将返回该岸的距离和所在国家/地区。

我正在使用来自以下国家/地区和海岸线的形状文件：http：//www.naturalearthdata.com/

一些考虑因素是不可访问的海洋极点为2688公里。因此，这将是距海岸的最大可能距离，可用于创建某种边界框。我想计算地球曲率（不是欧几里得）的解释，例如Haversine或Vincenty方法。

为此，我开始查看scipy.spatial.cKDTree，但这不适用于Haversine距离度量。另一方面，sklearn.neighbors.BallTree确实允许使用Haversine距离度量标准，但我无法使其正常工作。这是我到目前为止的代码。注意，理想情况下，应将功能向量化。

############################# 解决方案 ################## #############

感谢您的所有输入，这就是我用Python解决的方法，包括下载相关形状文件的功能，需要一些清洁

import os
import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

import shapely as sp
import cartopy.io.shapereader as shpreader
import ssl
import urllib.request
import zipfile

from shutil import rmtree
from dbfread import DBF
from scipy import spatial
from sklearn.neighbors import NearestNeighbors, BallTree
from pyproj import Proj, transform

from math import *

coastline = np.load(os.path.join(os.path.dirname(__file__),
                    '../data/shape_files/coast_coords_10m.npy'))

ports = np.load(os.path.join(os.path.dirname(__file__),
                '../data/shape_files/ports_coords.npy'))

def extract_geom_meta(country):
    '''
    extract from each geometry the name of the country
    and the geom_point data. The output will be a list
    of tuples and the country name as the last element.
    '''
    geoms = country.geometry
    coords = np.empty(shape=[0, 2])
    for geom in geoms:
        coords = np.append(coords, geom.exterior.coords, axis = 0)

    country_name = country.attributes["ADMIN"]
    return [coords, country_name]

def save_coastline_shape_file():
    '''
    store shp files locally, this functions will download
    shapefiles for the whole planet.
    '''
    ne_earth = shpreader.natural_earth(resolution = '10m',
                                       category = 'cultural',
                                       name='admin_0_countries')
    reader = shpreader.Reader(ne_earth)
    countries = reader.records()
    # extract and create separate objects
    world_geoms = [extract_geom_meta(country) for country in countries]
    coords_countries = np.vstack([[np.array(x[:-1]), x[-1]]
                                    for x in world_geoms])
    coastline = np.save(os.path.join(os.path.dirname(__file__),
                        '../data/shape_files/coast_coords_10m.npy')
                        , coords_countries)
    print('Saving coordinates (...)')

def distance_to_shore(lon, lat):
    '''
    This function will create a numpy array of distances
    to shore. It will contain and ID for AIS points and
    the distance to the nearest coastline point.
    '''
    coastline_coords = np.vstack([np.flip(x[0][0], axis=1) for x in coastline])
    countries = np.hstack([np.repeat(str(x[1]), len(x[0][0])) for x in coastline])
    tree = BallTree(np.radians(coastline_coords), metric='haversine')
    coords = pd.concat([np.radians(lat), np.radians(lon)], axis=1)
    dist, ind = tree.query(coords, k=1)
    df_distance_to_shore = pd.Series(dist.flatten()*6371, name='distance_to_shore')
    df_countries = pd.Series(countries[ind].flatten(), name='shore_country')
    return pd.concat([df_distance_to_shore, df_countries], axis=1)

Answer 1

解决此问题的有效方法是将测地距离作为度量标准（将度量标准满足三角形不等式很重要）， 将所有海岸点存储到有利的位置树中。然后，对于每个船只，您可以查询VP树以找到闭合点。

如果有M个海岸点和N个船只。然后，构造VP树的时间需要M log M距离计算。每个查询都需要对数M距离的计算。椭球的距离计算大约需要2.5μs。因此，总时间为（M + N）log M ×2.5μs。

这是使用我的库GeographicLib（版本1.47或更高版本）执行此计算的代码。这只是为NearestNeighbor类给出的示例的精简版本。

// Example of using the GeographicLib::NearestNeighbor class.  Read lon/lat
// points for coast from coast.txt and lon/lat for vessels from vessels.txt.
// For each vessel, print to standard output: the index for the closest point
// on coast and the distance to it.

// This requires GeographicLib version 1.47 or later.

// Compile/link with, e.g.,
// g++ -I/usr/local/include -lGeographic -L/usr/local/bin -Wl,-rpath=/usr/local/lib -o coast coast.cpp

// Run time for 30000 coast points and 46217 vessels is 3 secs.

#include <iostream>
#include <exception>
#include <vector>
#include <fstream>

#include <GeographicLib/NearestNeighbor.hpp>
#include <GeographicLib/Geodesic.hpp>

using namespace std;
using namespace GeographicLib;

// A structure to hold a geographic coordinate.
struct pos {
  double _lat, _lon;
  pos(double lat = 0, double lon = 0) : _lat(lat), _lon(lon) {}
};

// A class to compute the distance between 2 positions.
class DistanceCalculator {
private:
  Geodesic _geod;
public:
  explicit DistanceCalculator(const Geodesic& geod) : _geod(geod) {}
  double operator() (const pos& a, const pos& b) const {
    double d;
    _geod.Inverse(a._lat, a._lon, b._lat, b._lon, d);
    if ( !(d >= 0) )
      // Catch illegal positions which result in d = NaN
      throw GeographicErr("distance doesn't satisfy d >= 0");
    return d;
  }
};

int main() {
  try {
    // Read in coast
    vector<pos> coast;
    double lat, lon;
    {
      ifstream is("coast.txt");
      if (!is.good())
        throw GeographicErr("coast.txt not readable");
      while (is >> lon >> lat)
        coast.push_back(pos(lat, lon));
      if (coast.size() == 0)
        throw GeographicErr("need at least one location");
    }

    // Define a distance function object
    DistanceCalculator distance(Geodesic::WGS84());

    // Create NearestNeighbor object
    NearestNeighbor<double, pos, DistanceCalculator>
      coastset(coast, distance);

    ifstream is("vessels.txt");
    double d;
    int count = 0;
    vector<int> k;
    while (is >> lon >> lat) {
      ++count;
      d = coastset.Search(coast, distance, pos(lat, lon), k);
      if (k.size() != 1)
          throw GeographicErr("unexpected number of results");
      cout << k[0] << " " << d << "\n";
    }
  }
  catch (const exception& e) {
    cerr << "Caught exception: " << e.what() << "\n";
    return 1;
  }
}

此示例在C ++中。要使用python，您需要找到VP树的python实现，然后可以使用 python版本的GeographicLib进行距离计算。

PS GeographicLib对满足三角形不等式的测地距离使用精确的算法。Vincenty方法无法收敛几乎对映点，因此不能满足三角形不等式。

附：这是python实现：安装vptree和geogelib

pip install vptree geographiclib

海岸点（lon，lat）在Coast.txt中；容器位置（lon，lat）位于vessel.txt中。跑

import numpy
import vptree
from geographiclib.geodesic import Geodesic

def geoddist(p1, p2):
  # p1 = [lon1, lat1] in degrees
  # p2 = [lon2, lat2] in degrees
  return Geodesic.WGS84.Inverse(p1[1], p1[0], p2[1], p2[0])['s12']

coast = vptree.VPTree(numpy.loadtxt('coast.txt'), geoddist, 8)
print('vessel closest-coast dist')
for v in numpy.loadtxt('vessels.txt'):
  c = coast.get_nearest_neighbor(v)
  print(list(v), list(c[1]), c[0])

对于30000个海岸点和46217艘船，这需要18分3秒。这比我预期的要长。构造树的时间为1分16秒。因此，总时间应为3分钟左右。

对于30000个海岸点和46217艘船，这需要4分钟（使用vptree 1.1.1版）。为了进行比较，使用GeographicLib C ++库的时间为3秒。

后来：我调查了为什么python vptree速度慢。设置树的距离计算次数对于GeographicLib的C ++实现和python vptree包是相同的：387248，大约M log M，对于M =30000。（此处的日志为2，我将存储桶大小设置为1这两种实现方式都可以简化比较。）对于C ++实现，每个容器查找的距离计算的平均数为14.7，接近预期值，log M = 14.9。但是，python实现的等效统计数据为108.9，这是7.4的一个大因素。

各种因素影响VP树的效率：华帝点的选择，搜索是如何排序的，这些因素为GeographicLib实现给定等的讨论在这里。我将对此ping python软件包的作者。

稍后：我已经提交了pull请求，该请求解决了python软件包vptree效率方面的主要问题。我的测试的CPU时间现在约为4分钟。每个查询的距离计算数量为16.7（接近于GeographicLib :: NearestNeighbor的数字，为14.7）。

Answer 2

这里的关键是您需要使用“大圆”（正交）距离计算，该计算旨在查找球体表面上两点之间的距离。虽然地球不是一个完美的球体，但这样的计算将使您非常接近（在 0.5% 以内），如果不够接近，可以应用非球形调整。

互联网上有很多关于这个公式的文档。您将需要寻找涉及 XYZ 而不是极坐标的封闭形式解决方案，或者将 GPS 坐标转换为极坐标（两者之一）。