Skip to content

Koordinaten umrechnen in R

Geocaching Vor knapp 10 Jahren (also in Computerzeitrechnung vor einer Ewigkeit) habe ich hier Excel-Tabellen gezeigt um Koordinaten zwischen verschiedenen Bezugssystemen umrechnen zu können.

In der letzten Zeit musste ich mich viel mit R beschäftigen. R erlebt derzeit eine immense Verbreitung, es ist eine einfache Programmiersprache die viel für Statistik verwendet wird. Man kann damit beachtliche Datenmengen umrechnen und analysieren. Eine große Sammlung an Packages mit denen die Sprache für den jeweiligen Zweck erweitert wird ermöglicht unter anderem auch das Rechnen mit Koordinaten.

Wie schön und einfach das geht kann ich hier an einem Beispiel aus der grünen Hölle zeigen. Zunächst einmal verwenden wir die Bibliotheken "sp" für räumliche Objekte und "geospehere" für geodätische Hauptaufgaben auf dem Ellipsoid.

library("sp")
library("geosphere")


Anschließend machen wir zwei Vektoren für die beiden Koordinaten, die in unterschiedlichen UTM-Zonen vorliegen, und erzeugen sogenannte SpatialPoints daraus:

a=cbind(710482,5654109)
ap=SpatialPoints(a, CRS("+proj=utm +zone=32"))
b=cbind(289511,5654098)
bp=SpatialPoints(b, CRS("+proj=utm +zone=33"))


Hier könnten es auch ganz andere Koordinaten sein, also z.B. Gauß-Krüger oder Soldner, das muss dann nur bei der Erzeugung der SpatialPoints angegeben werden. Die Umrechnung in ein gemeinsames Bezugssystem ist nun äußerst einfach und mit je einem Befehl erledigt:

ae=spTransform(ap, CRS("+proj=longlat +ellps=WGS84"))
be=spTransform(bp, CRS("+proj=longlat +ellps=WGS84"))


Dabei wäre auch die Angabe der beliebten EPSG-Codes möglich (also in diesem Fall EPSG:4326). Für alle die sich mit Geocaching beschäftigen: Das sind die Koordinaten so wie sie dort verwendet werden, also Länger und Breite auf dem Ellipsoid.

Danach kann man sehr einfach mit diesem Koordinaten rechnen, also z.B. den Abstand und den Richtungswinkel:

d=distVincentyEllipsoid(ae,be)
print(d)
b=bearing(ae,be)
print(b)


Auch die bei Geocachern beliebte Waypoint-Projektion (die eher ein polares Anhängen ist) kann man so einfach durchführen.

Fastfood an der Autobahn fürs Garmin

OpenStreetMap So langsam wirds draußen wieder wärmer und es gibt die eine oder andere Strecke auf der Autobahn zurückzulegen. Schön auch dass die großen Fastfood-Ketten dort auch ihren Platz gefunden haben, mit bekannten Produkten zu zivilen Preisen. Gerade beim Geocaching hat die "Goldene Möve" schon fast Tradition.

Nur die Suche nach dem nächsten Fastfood ist gar nicht so einfach: Das Navi kennt natürlich alles, aber schlägt sortiert nach Entfernung vor. Ich suche in diesen Fällen aber in der Nähe einer der nächsten Autobahn-Abfahrten, also z.B. einen Autohof der nicht so weit ab von der Autobahn-Abfahrt liegt.

Hier war also eine neue POI-Datei fürs Garmin gefragt. Wer nur die POI haben will, kann hier die Datei motorway_fastfood.gpi herunterladen und einfach in den POI-Ordner kopieren.

Um die Fastfood-POI an der Autobahn aus der OpenStreetMap abzuleiten, verwende ich in diesem Beispiel einen deutschen Auszug der OSM, der mit osm2pgsql in eine PostGIS-Datenbank eingespielt wurde.

Alle Fastfood-Läden sind einfach selektiert:

select osm_id, way as coords, name from planet_osm_point where amenity='fast_food'
union
select osm_id, st_centroid(way) as coords, name from planet_osm_polygon where amenity='fast_food';


Auch die Autobahnen liegen hier bereits vor:

select ref, way from planet_osm_roads where highway='motorway'


Beides wird mit einem "spatial join" den PostGIS anbietet verknotet. Man kann mit dem "on" nicht nur Schlüssel matchen wie sonst in SQL üblich, sondern auch einen maximalen Abstand abfragen. Komplett sieht die Abfrage dann so aus:
copy (
select distinct on (osm_id)
st_x(st_transform(poi.coords,4326)) as lon, st_y(st_transform(poi.coords,4326)) as lat, concat(bab.ref,' ',poi.name) as name
from (
select osm_id, way as coords, name from planet_osm_point where amenity='fast_food'
union
select osm_id, st_centroid(way) as coords, name from planet_osm_polygon where amenity='fast_food'
) as poi
inner join (
select ref, way from planet_osm_roads where highway='motorway'
) as bab
on st_distance(bab.way,poi.coords)<400 order by osm_id
) to '/tmp/motorway_fastfood.csv' CSV header;


Damit werden die Daten auch gleich in eine CSV-Datei kopiert, die dann noch mit GPS-Babel auf der Shell in ein Garmin-taugliches GPI-File mit einem Burger als Signatur konvertiert werden kann:
wget http://openclipart.org/image/800px/svg_to_png/9075/Gerald_G_Fast_Food_Lunch_Dinner_%28FF_Menu%29_1.png -O motorway_fastfood.png

convert motorway_fastfood.png -compress none -type truecolor -resize 18x18 motorway_fastfood.bmp

gpsbabel -c utf8 -i unicsv -f "motorway_fastfood.csv" -o garmin_gpi,category="fastfood",bitmap="motorway_fastfood.bmp",position -F "motorway_fastfood.gpi"


Einen Haken hat die Sache noch: In großen Städte sind an den Abfahrten in der Stadt sehr viele Fastfood-POI die natürlich eigentlich nicht die gesuchten sind, weil sie sich nicht für die kurze Pause eignen. Für mich ist das kein Problem weil ich diese POI-Sammlung ohnehin nur in weniger dicht besiedelten Gebieten benötige.

Update: Ein erster Praxistest ergab, dass die SQL-Abfrage mit den Ausfahrten nicht so günstig war und daher wird nun rund um die Autobahn an sich gesucht. Sehr schön dass man im Nüvi damit jetzt auch nach dem Namen der Autobahn suchen kann auf der man gerade fährt!

Die Erstellung dieser POI-Datei ist gleichzeitig auch ein schönes Beispiel wie man allgemein POI aus der OpenStreetMap aufs Garmin bekommen kann. Lasst mich in den Kommentaren wissen wie Ihr mit der Datei zurecht kommt!

Frosch-Post aus Seattle

Geocaching Da freut sich aber der Ingress-Frosch: Groundspeak bedankt sich für meinen sportlichen Einsatz bei der Wahl zum Geocacher des Monats mit einem kleinen Geschenk-Package.

Der Brief aus Seattle war heute mittag in meiner Mailbox im Hauseingang. Mit dabei war Kermit, ähh, Signal der Frosch, in Form der 2012er Lackey Geocoin, der sofort der Bosch'schen Perforationskunst unterzogen wurde. So kann ich meine neue Frosch-Geocoin an die Kette (oder besser den Lanyard) legen.

Es wurde auch langsam Zeit, meinen 17er Schlüssel den ich seit einigen Jahren auf Events ausführe durch einen Frischen zu ersetzen.

GST049

Passend zu meiner neuen Freizeitbeschäftigung ein grüner Frosch. Perfekt. Der wird heute abend auf Mic@'s Seifenblasen-Event sofort eingeweiht.
tweetbackcheck