Koordinaten umwandeln von Dezimalgrad in projizierte Koordinaten

Wie kann ich Koordinaten in Dezimalgrad, z.B. 52.5185584/13.3761286 (Breite/Länge), in das umwandeln, was JOSM als projizierte Koordinaten bezeichnet, also 1.0809483/0.2334575 für das Beispiel? Das ganze soll in einem Script erfolgen. Wahrscheinlich fehlen mir nur die richtigen Begriffe, denn mit den mir bekannten komme ich irgendwie nicht weiter. cs2cs aus der proj-Bibliothek (http://trac.osgeo.org/proj/) scheint so was zu können, aber was sind da die richtigen Optionen? Hintergrund ist, dass ich aus den Koordinaten eines Luftbildes eine PicLayer-Kalibrierungsdatei generieren will.

Das dürfte nicht möglich sein. Die Erde ist eine Kugel. Dein Monitor ist keine Kugel. Es gibt 100erte Projektionsmethoden wie man jetzt die Karte anzeigen kann, die optimale ist “Merkator”.

Schaue z.B. mal auf einem Globus vom Äquator auf Deutschland. Die Stirn muss am Äquator aufliegen. Jetzt ist Deutschland viel breiter als hoch, weil Du schief schaust. Die Projektionsmethode Merkator schaut z.B. immer “von Oben” damit dieses eine Problem nicht vorhanden ist.

JOSM speichert jetzt zu jedem Node die XY-Koordinaten wie es am Monitor aus sieht. Möchtest Du mit einem Plugin ein Rechteck machen, sucht sich das Plugin die XY-Koordinaten und macht sie rechteckig. Das Problem: Stellst Du auf eine andere Projektion um, ist es plötzlich nicht mehr rechteckig. Damit die Plugins es einfacher haben errechnet JOSM aber diese Werte.

Wenn Du diese Werte meinst: Die kannst Du nicht errechnen, da sie Abhängig von der JOSM-Ebene sind und bei jedem JOSM-Neustart anders sind.

Sehr schöne Erklärung!
Würde mich freuen wenn die im Wiki auch auftaucht. (z.B. http://wiki.openstreetmap.org/wiki/Projektionen))
Am Anfang tat ich mich auch richtig schwer das zu verstehen.

Gruß

Vielleicht hilft das weiter Slippy map tilenames: Formeln, Programmcode in verschiedenen Sprachen.

DIE optimale Projektion gibt es nicht. Es gibt immer irgendwo einen Haken. Bei Merkator ist es, daß sie nicht flächengetreu ist. Objekte nahe den Polen werden größer angezeigt, als am Äquator. Deswegen sieht Grönland auch so groß aus wie ganz Afrika. Wenn man das mal auf einem Globus nachsieht, dann stellt man fest, daß Grönland viel viel kleiner ist.

Christian

@brogo: Ja, das stimmt. Ich ergänze meinen Satz zu “… die optimale Projektionsmethode für Deutschland und Openstreetmap ist für Standardanwendungen Merkator”

@godofglow: Fühl Dich frei, den Text zu verwenden. Ich erhebe keine Urheberrechte drauf und erdenke auch nicht, jemanden zu verklagen :wink:

Dass die Erde eine Kugel ist, ist etwas sehr vereinfacht. Ellipsoid passt besser. Geoid wäre noch genauer. Die Auswirkungen der verschiedenen Projektionen kann man auch schön im JOSM sehen, wenn slippymap aktiv ist und man gaaanz weit rauszoomt.

Mir geht es aber mehr um Koordinaten. Wenn ich einen GPX-Track im JOSM lade, dann werden die einzelnen Wegpunkte doch immer an der Stelle angezeigt, zu der sie gehören, unabhängig von der Projektion. Klar, die Abstände der Punkte auf dem Monitor verändern sich mit der Projektionsmethode und die Wege werden mehr oder weniger verzerrt, aber es wird deswegen nichts falsch angezeigt. Im GPX steht jeder Punkt als lat/lon-Wertepaar, z.B. . Diese Zahlen werden im JOSM auch als Breit und Länge angezeigt, wenn ich unter Karten-Einstellungen Merkator als Projektionsmethode und “Koordinaten anzeigen als” Dezimalgrad einstelle. Wenn ich “Koordinaten anzeigen als” “Projizierte Koordinaten” einstelle, dann werden diese Koordinaten doch sicher auch eindeutig sein. Allerdings ändern sich die Werte der “projizierten Koordinaten” tatsächlich mit der Projektion. Und ich suche nun nach einem Weg die Latitude/Longitude Werte des GPX, bzw. das, was bei Projektion Merkator Dezimalgrad heißt, in das Format umzuwandeln, was bei Projektion Merkator als “Projizierte Koordinaten” im JOSM bezeichnet wird. Ich muss wahrscheinlich mal in der englischen Version nachsehen, um die richtigen Begriffe zu finden…

Die Seite mit den Tilenames hat nicht weiter geholfen, da dort immer nur mit lon/lat hantiert wird.

Ich konnte erst nicht viel damit anfangen, aber letztendlich war das genau der richtige Tip. Die Transformation geht mit Perl so:


use Math::Trig;
$lonproj = $lon / 180.0 * pi;
$latproj = log(tan(deg2rad($lat)) + sec(deg2rad($lat)));

Ist allerdings nur für positive Koordinaten verifiziert. Vielen Dank für die Hilfe. Die damit generierte PicLayer-Kalibrierungsdatei sieht auch ganz gut aus. Die Koordinate, die man dort angeben muss, ist genau der Bildmittelpunkt. An der genauen Bedeutung der Scale knobel ich noch.

Das hört sich ja sehr interessant an! Wäre es möchlich, diese Skript (oder was es auch ist) auch zu bekommen? Wäre echt super!

Ist möglich. Die Frage ist nur, was du damit machen willst bzw. was der Input des Skriptes sein soll. Bei mir geht es um einen Haufen kleiner Kacheln, die ich zu mehreren größeren Kacheln zusammensetzen und dann in JOSM laden will. Also wahrscheinlich eher was spezielles. Ist auch noch nicht fertig.

Als “Input” für das Skript stehen (georeferenzierte) “Landkartenskizzen” eines Freundes zu Verfügung ( er hat der Verwendung für OSM zugestimmt) und diese sind erstaunlich gut, zumindest mit bisl Ortskenntnis. Ich habe sie als einfache A4 Blätter bekommen. Das Scannen war noch ein klaks, dass georeferenzieren dann schon um längen schwerer und nervenraubend und jetzt will ich es nicht nochmal “per Hand” für JOSM machen. Also wenn du das Skript “fertig” hast (oder auch jetzt schon) würd ich mich über die Bereitstellung freuen.

Hier kommt das Skript. Es hat allerdings den Schwachpunkt, dass die Scale nicht ganz stimmt. Ich habe schon alles probiert, was mir eingefallen ist, aber die Bildgröße ist nicht ganz korrekt (zu groß). Den Teil mit Koordinaten und Bildauflösung bestimmen habe ich weggelassen, da der Weg je nach Quelle variiert. In meinem Fall steht das alles schon im Dateinamen.


#! /usr/bin/env perl
###
### makecal.pl - Generate a PicLayer calibration file based on two coordinates.

### Input is two coordinates of opposite image corners and the image
### resolution.
### Lower left corner:
$latx = 52.5180789;
$lonx = 13.3753334;
### Upper right corner:
$laty = 52.5190411;
$lony = 13.3769147;
### Width (left to tight) in pixels:
$width = 500;
### Height (top to bottom) in pixels:
$height = 500;
### Name if the image file:
$imgFile = "image.jpg";

use POSIX;
use Math::Trig;
use warnings;

### Use the 'geod' tool from the PROJ.4 package to calculate the
### distance between two coordinates?  Some almost correct math will
### be used otherwise.
$useGeod = 1;

my ($xmin, $ymin) = latlon2xy($latx, $lonx);
my ($xmax, $ymax) = latlon2xy($laty, $lony);
my $xpos = ($xmin + $xmax) / 2.0;
my $ypos = ($ymin + $ymax) / 2.0;
my $calScale = getScaleWH($latx, $lonx, $laty, $lony, $width, $height);
open(my $calFd, ">", $imgFile.".cal") or die $!;
print $calFd <<EOT;
INITIAL_SCALE=$calScale
SCALEX=1.0
SCALEY=1.0
ANGLE=0.0
INITIAL_POS_X=$xpos
INITIAL_POS_y=$ypos
POSITION_X=$xpos
POSITION_Y=$ypos
EOT
close($calFd);

######################################################

### Convert latitude and longitude into x/y coordinates.
sub latlon2xy {
  my($lat, $lon) = @_;
  my $lonproj = deg2rad($lon);
  my $latproj = log(tan(deg2rad($lat)) + sec(deg2rad($lat)));
  return($lonproj, $latproj);
}

### Calculate the distance between two coordinates in meters.
sub getDist {
  my ($latx, $lonx, $laty, $lony) = @_;
  my $dist;
  ### Source: http://www.kompf.de/gps/distcalc.html
  if ($useGeod) {
    ### WGS84 is the same as EPSG:4326.
    open(my $geod, "echo $latx $lonx $laty $lony"
     ." | geod +ellps=WGS84 +units=m -I |") or die $!;
    my $data = <$geod>;
    close($geod);
    $data =~ m/\S+\s+\S+\s+(\S+)/;
    $dist = $1;
  }
  else {
    $dist = 6378137.0 * acos(sin(deg2rad($latx)) * sin(deg2rad($laty))
                 + (cos(deg2rad($latx)) * cos(deg2rad($laty))
                * cos(deg2rad($lony - $lonx))));
  }
  return($dist);
}

### Calculate the scale based on the coordinates of the lower left and
### upper right corner of an image and the width and hight of that
### image.
sub getScaleWH {
  my ($latx, $lonx, $laty, $lony, $width, $height) = @_;
  my $dist = getDist($latx, $lonx, $laty, $lony);
  my $pixDist = sqrt($width**2 + $height**2);
  my $scale = $dist / $pixDist * 100.0;
  return($scale);
}

Hi Holger,

mmh, das Aufrufen von geod ist nicht ganz so schön :wink:
Schon mal
use Geo::Ellipsoid;
probiert?

Ciao,
Frank

Nein, kannte ich noch nicht. Ist leider bei meiner Disti (OpenSuse) nicht dabei. Muss ich erstmal nachsehen, wie das mit dem manuellen Installieren von Modulen als nicht-root geht. Es gibt allerdings Geo::Proj4 als Modul. Wäre das auch schöner?

$ cpan
cpan> o conf makepl_arg “INSTALL_BASE=/home/myuser”
cpan> o conf commit
cpan> install Geo::Ellipsoid


use lib ‘/home/myuser/lib/perl5’;

oder so ähnlich :wink:

Geo::Proj4 hat das überhaupt eine Routine zur Berechnung des Abstandes von Punkten auf Ellipsoid-Fächen?

Ciao,
Frank

Vor allem: Hat da schon mal jemand unter Windows zum Laufen gebracht?

Gruß,
ajoessen

Das mit geod ist alles Käse. Einfach $useGeod = 0 setzen und die Berechnung stimmt. Zumindest in meinem Fall passt es dann sehr gut. Ich hatte versucht möglichst genau zu rechnen, aber das ist hier verkehrt.

Geo::Proj4 scheint nur cs2cs zu implementieren, das geht also nicht.

@ajoessen: Du kannst einfach die Haversine-Formel nehmen, da diese vermutlich nur um wenige hundert Meter abweicht. Wenn Du es aber genauer haben willst, such’ einfach nach der Vincenty-Implementierung. Ich habe eine in Python und eine in PHP gemacht - Perl kann ich aber leider nicht. Das Netz ist aber voll damit.

Ich auch nicht. Ich möchte halt nur das mapgen.pl unter Windows zum Laufen bringen.
Der Entwickler von Goe::Porj4 ist leider Windows-Verweigerer.
Ich hab meinen Kenntnisstand auf:
http://wiki.openstreetmap.org/wiki/Talk:Mapgen.pl#Geo::Proj4_on_Windows
und hoffe auf Besserung.

Für meine eigene Programmierung hab ich diverse Formeln, die es in allen Sprachen tun :wink:

Gruß,
ajoessen

Hi,

naja, recht viel technische Details sind dort nicht zu finden. Hatte es mal mit Strawberry Perl probiert,
das hat ja einen C-compiler dabei.
Kompiliert wird’s auch noch, leider steigt dann der linker aus :frowning:

Ciao,
Frank