Berechnen Sie den Abstand zu Ufer oder Küstenlinie für ein Schiff

Für einen Datensatz von 200M GPS (lon, lat) Koordinaten von Schiffen möchte ich eine ungefähre Distanz zum nächstgelegenen Land oder Küstenlinie berechnen, als eine Funktion namens distance_to_shore, die die Distanz und das Land dieses Ufers zurückgeben wird.

Ich verwende eine Formdatei von Ländergrenzen und Küsten von: http://www.naturalearthdata.com/

Einige Überlegungen sind, dass der ozeanische Pol der Unzugänglichkeit 2688 km ist. So wäre dies die maximal mögliche Entfernung vom Ufer, das könnte verwendet werden, um irgendeine Art von Bounding Box zu schaffen. Ich möchte die Buchführung für die Krümmung der Erde (nicht Euklidisch), zB Haversine oder Vincenty Methode, berechnen.

Dafür habe ich angefangen, scipy.spatial.cKDTree zu betrachten, aber das erlaubt nicht Haversine Distanz metrisch. Auf der anderen Seite die Sklearn.neighbors.BallTree, erlaubt Haversine Abstand metrisch, aber ich kann es nicht zu arbeiten. Hier ist der Code, den ich bisher habe. NB sollte die Funktion idealerweise vektorisiert werden.

def proj_arr(points, proj_to): """ Project geographic co-ordinates to get cartesian x,y Transform(origin|destination|lon|lat) to meters. """ inproj = Proj(init='epsg:4326') outproj = Proj(init=proj_to) func = lambda x: transform(inproj, outproj, x[0], x[1]) return np.array(list(map(func, points))) def distance_to_shore(lon, lat, country_name=False): ''' 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. ''' #lon = pd.Series(-79.049683) #lat = pd.Series(4.328306) coastline_coords = np.vstack([x[0] for x in coastline]) countries = np.hstack([np.repeat(str(y[1]), len(y[0])) for y in coastline]) coords = pd.concat([lon, lat], axis=1) coords_matrix = coords.as_matrix([coords.columns[0:2]]) #Project to meters using 'proj_arr' function and calculate distance coast_proj = proj_arr(coastline_coords, 'epsg:3410') coords_proj = proj_arr(coords_matrix, 'epsg:3410') #tree = BallTree(coastline_coords, metric='haversine') #tree = BallTree(coastline_coords, metric='haversine') #tree= cKDTree(coast_proj) #distance, index = tree.query(coords) #distance_upper_bound= 2688 km distance, index = spatial.cKDTree(coast_proj).query(coords_proj) sea = coords_matrix land = coastline_coords[index] distance = great_circle(sea,land).kilometers df_distance_to_shore = pd.Series(distance, name='distance_to_shore') df_countries = pd.Series(countries[index], name='shore_country') return pd.concat([df_distance_to_shore, df_countries], axis=1) 

Einige dieser Zeilen in der distance_to_shore-Funktion sind überflüssig, zeigen aber den anderen Ansatz, den ich versucht habe …

Vielen Dank

    3 Solutions collect form web for “Berechnen Sie den Abstand zu Ufer oder Küstenlinie für ein Schiff”

    Der effiziente Weg, dieses Problem zu lösen, besteht darin, alle Ihre Küstenpunkte in einen Aussichtspunktbaum zu bringen, indem man den geodätischen Abstand als Metrik verwendet (es ist wichtig, dass die Metrik die Dreiecksungleichheit erfüllt). Dann können Sie für jedes Schiff den VP-Baum abfragen, um den geschlossenen Punkt zu finden.

    Wenn es M Küstenpunkte und N Schiffe gibt. Dann benötigt die Zeit zum Aufbau des VP-Baums M log M Abstandsberechnungen. Jede Abfrage erfordert log M Abstand Berechnungen. Eine Abstandsberechnung für das Ellipsoid dauert ca. 2,5 μs. So ist die Gesamtzeit ( M + N ) log M × 2,5 μs.

    Hier ist Code mit meiner Bibliothek GeographicLib (Version 1.47 oder höher), um diese Berechnung durchzuführen. Dies ist nur eine abgespeckte Version des Beispiels für die NearestNeighbor-Klasse gegeben .

     // 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, eg, // 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; } } 

    Dieses Beispiel ist in C ++. Um Python verwenden zu können, musst du eine Python-Implementierung von VP-Bäumen finden und dann kannst du die Python-Version von GeographicLib für die Entfernungsberechnungen verwenden.

    PS GeographicLib verwendet einen genauen Algorithmus für den geodätischen Abstand, der die Dreiecksungleichung erfüllt. Die Vincenty-Methode fällt nicht für nahezu antipodale Punkte zusammen und erfüllt damit nicht die Dreiecksungleichung.

    ADDENDUM : Hier ist die Python-Implementierung: Installiere vptree und geographiclib

     pip install vptree geographiclib 

    Küstenpunkte (lon, lat) sind in coast.txt; Schiffspositionen (lon, lat) sind in Gefäßen.txt. Lauf

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

    Für 30000 Küstenpunkte und 46217 Schiffe dauert dies 18 min 3 sec. Das ist länger als ich erwartet hatte. Die Zeit, den Baum zu bauen, ist 1 min 16 Sek. Also die Gesamtzeit sollte ca. 3 min betragen.

    Zum Vergleich ist die Zeit mit der GeographicLib C ++ – Bibliothek 3 Sek.

    LATER : Ich sah hin, warum die Python Vptree ist langsam. Die Anzahl der Distanzberechnungen zum Einrichten des Baumes ist für die C ++ – Implementierung von GeographicLib und das python vptree-Paket gleich: 387248, was etwa M log M für M = 30000 ist. (Hier sind die Protokolle Basis 2 und ich setze die Eimergröße auf 1 für Beide Implementierungen, um Vergleiche zu erleichtern.) Die mittlere Anzahl der Entfernungsberechnungen für jede Schiffssuche nach der C ++ – Implementierung beträgt 14,7, die nahe dem erwarteten Wert liegt, log M = 14,9. Allerdings ist die äquivalente Statistik für die Python-Implementierung 108,9, ein Faktor für 7,4 größer.

    Verschiedene Faktoren beeinflussen die Effizienz des VP-Baumes: die Wahl der Aussichtspunkte, wie die Suche bestellt wird usw. Eine Diskussion dieser Überlegungen für die GeographicLib-Implementierung ist hier gegeben . Ich werde den Autor des Python-Pakets dazu picken.

    Der Schlüssel hier ist, dass Sie die "großen Kreis" (orthodrom) Abstand Berechnungen verwenden müssen, die entworfen sind, um den Abstand zwischen zwei Punkten auf der Oberfläche einer Kugel zu finden . Obwohl die Erde nicht eine perfekte Sphäre ist, werden solche Berechnungen Sie sehr nah (bis zu 0,5%) erhalten, und nicht-sphärische Anpassungen können angewendet werden, wenn dies nicht nahe genug ist.

    Es gibt viele Dokumentationen dieser Formel im Internet. Sie wollen nach geschlossenen Formular-Lösungen suchen, die XYZ anstelle von Polarkoordinaten beinhalten oder Ihre GPS-Koordinaten in Polar, eine der beiden umwandeln.

    Sie benötigen eine grosse Kreisberechnungsformel. Diese werden manchmal als Sphärisches Cosinusgesetz , Haversine oder Vincenty , Formeln bezeichnet.

    Sie können dann die Entfernung von jedem Schiff bis zum nächsten Punkt in Ihrem Küsten-Korpus berechnen. Es ist oft hilfreich, eine Bounding-Box-Berechnung zu verwenden, um irrelevante Punkte auszuschließen, bevor sie die ganze Great Circle-Formel auf sie laufen lassen.

    Wenn du dein Küstenkorpus konstruierst, musst du vielleicht die Interpolation verwenden, um zusätzliche Küstenspitzen hinzuzufügen, wenn deine Rohküstendaten lange Segmente haben. Das ist, weil Sie berechnen Entfernung zum nächsten Punkt nicht das nächste Segment . Schau auf große Kreisinterpolation.

    Wenn deine Schiffe in der Nähe von Pole sind (gut, in der Nähe des Nordpols, da man sieht, wie der Südpol auf dem Land ist) werden die Dinge mit den Standard-Großen Kreisformeln und den umgrenzenden Rechtecken funky werden. Sie sollten wahrscheinlich die Vincenty Formulierung in diesem Fall verwenden.

    Hier ist ein Aufschreiben über die Verwendung eines DBMS mit Indizierung für diese Art von Zweck. https://www.plumislandmedia.net/mysql/haversine-mysql-nearest-loc/

    Wenn Sie NOAA-Chart Level Genauigkeit benötigen, müssen Sie wahrscheinlich über Universal Transverse Mercator Projektionen zu lernen. Das ist über den Rahmen einer Stack Overflow Antwort hinaus.

    Python ist die beste Programmiersprache der Welt.