Kartographie (yukterez.net Backup)

Mollweide

Welt- oder Himmelskarten liegen gängigerweise entweder in der Plattkarten- oder der Mollweide-Projektion vor. Der folgende Artikel zeigt wie man einen Kartennetzentwurf in den anderen umwandelt, und wie man ihn auf eine Kugel projiziert. Das in den Beispielen verwendete Rohmaterial ist das Wikipedia-Bild Equirectangular-projection.jpg

Zuerst wird die Funktion für die Koordinatentransformation definiert und ihr der Name Eq2Mol zugewiesen:

Code:

Eq2Mol[{x_, y_}] := {Lmol[x, y], Bmol[y]};
θ[y_, rad_: 1] := ArcSin[y];
Bmol[y_] := ArcSin[(2 θ[y] + Sin[2 θ[y]])/π];
Lmol[x_, y_] :=  0 + π x/(2 Cos[θ[y]]);
(* Equirectangular <> Mollweide *)

Anstatt der 0 kann auch der Winkel um den der Längengrad verschoben werden soll in Radianten eingegeben werden. Dann wird eine Karte geladen, in diesem Fall eine Plattkarte, und dem Bild eine beliebige Variable, hier Eq für Equirectangular, zugewiesen:

Code:

Eq = Import["http://upload.wikimedia.org/wikipedia/commons/e/ea/Equirectangular-projection.jpg"]
(* Plattkarte laden *)

Die Plattkarte Eq wird ins Mollweide-Format umgewandelt und das generierte Bild als Mw gespeichert:

Code:

Mw = ImageTransformation[Eq, Eq2Mol,
  DataRange -> {{-π, π}, {-π/2, π/2}},
  PlotRange -> {{-2, 2}, {-1, 1}}]
  (* Equirectangular -> Mollweide *)

Und umgekehrt von Mollweide wieder zurück nach Equirectangular als eQ:

Code:

eQ = ImageForwardTransformation[Mw, Eq2Mol,
  DataRange -> {{-2, 2}, {-1, 1}},
  PlotRange -> All]
  (* Mollweide -> Equirectangular *)

Da es sich um eine Pixeloperation handelt werden die Qualitätsverluste desto höher je öfter man von einer Projektion in die andere und wieder zurück transformiert.

Hammer-Aitoff

Die Funktion um von der Equirectangular- in die Hammer-Aitoff-Projektion zu transformieren lautet:

Code:

Eq2HA[{x_, y_}] := {Lha[x, y], Bha[x, y]};
ζ[x_, y_] := Sqrt[1 - x^2/16 - y^2/4];
Bha[x_, y_] := ArcSin[ζ[x, y] y];
Lha[x_, y_] := 2 ArcTan[ζ[x, y] x/(2 (2 ζ[x, y]^2 - 1))];
(* Equirectangular >> Hammer Aitoff *)

Das verwendete Rohmaterial ist wieder dieses:

Code:

Eq = Import["C:\\Users\\Yukterez\\Desktop\\Equirectangular-projection.jpg"]
(* Plattkarte laden *)

Die Plattkarte Eq wird nun in die Hammer-Aitoff Projektion Ha umgewandelt. Die doppelt projizierten Bereiche an den Rändern können bei Bedarf mit einer schwarzen oder weissen Ellipsenumhüllung abgedeckt werden:

Code:

Ha = ImageTransformation[Eq, Eq2HA,
  DataRange -> {{-π, π}, {-π/2, π/2}},
  PlotRange -> {{-2 Sqrt[2], 2 Sqrt[2]}, {-Sqrt[2], Sqrt[2]}}]
  (* Equirectangular -> Hammer Aitoff *)

Für die Rückverwandlung ins Plattkartenformat wird die Transformation

Code:

HA2Eq[{x_, y_}] := {X[x, y], Y[x, y]};
X[x_, y_] := (2 Sqrt[2] Cos[y] Sin[x/2])/Sqrt[1 + Cos[y] Cos[x/2]];
Y[x_, y_] := (Sqrt[2] Sin[y])/Sqrt[1 + Cos[y] Cos[x/2]];
(* Hammer Aitoff >> Equirectangular *)

benötigt (da bei der Transformation in das Hammer-Aitoff-Format die selben Koordinaten sowohl ins Bild als auch an den Rand projiziert werden kann bei einer Invertierung der Funktion nicht wie sonst üblich die ImageTransformation[] in eine ImageForwardTransformation[] überführt werden).

Code:

eq = ImageTransformation[Ha, HA2Eq,
 DataRange -> {{-2 Sqrt[2], 2 Sqrt[2]}, {-Sqrt[2], Sqrt[2]}},
 PlotRange -> {{-π, π}, {-π/2, π/2}}]
 (* Hammer Aitoff -> Equirectangular *)

Quincunx

Um eine Plattkarte ins Quincunx-Format zu transformieren wird wie immer zuerst die Plattkarte geladen:

Code:

Eq = Import["C:\\Users\\Yukterez\\Desktop\\Equirectangular-projection.jpg"]
(* Plattkarte laden *)

Danach wird die Transformationsregel implementiert:

Code:

Eq2Pq[{x_, y_}] := {Arg[JacobiCN[x + I y, 1/2]], 2 ArcCot[Abs[JacobiCN[x + I y, 1/2]]]};
R = EllipticK[1/2];
(* Quincunx <> Equirectangular *)

oder wenn es um reine Pixeltransformationen geht alternativ dazu:

Code:

j[ε_] := Max[0, Ceiling[Log2[4 Abs[ε]]]];
ξ[ε_] := (ε 2^-j[ε])^2;
fPq[ε_] := Nest[-(((#^2+2) #^2-1)/((#^2-2) #^2-1)) &, (1-ξ[ε]/4 (1+ξ[ε]/30 (1+ξ[ε]/8)))/(1 + ξ[ε]/4 (1 - ξ[ε]/30 (1 - ξ[ε]/8))), j[ε]];
Eq2Pq[{x_, y_}] := {Arg[fPq[x + I y]], 2 ArcCot[Abs[fPq[x + I y]]]};
R = EllipticK[1/2];
(* Quincunx <> Equirectangular *)

Und angewandt:

Code:

Pq = ImageTransformation[Eq, Eq2Pq,
  DataRange -> {{-π, +π}, {-0, +π}}, PlotRange -> {{0, 2 R}, {-R, R}}]
  (* Equirectangular -> Quincunx *)

Da das so gezeichnete Muster ein wiederholendes ist kann es endlos fortgesetzt werden:

Code:

pq = ImageTransformation[Eq, Eq2Pq,
  DataRange -> {{-π, +π}, {-0, +π}}, PlotRange -> {{0, 4 R}, {-R, 3 R}}]
  (* Equirectangular -> Quincunx *)

Zurücktransformiert wird mit

Code:

eP = ImageForwardTransformation[pq, Eq2Pq,
  DataRange -> {{0, 4 R}, {-R, 3 R}}, PlotRange -> {{-π, +π}, {-0, +π}}]
  (* Quincunx -> Equirectangular *)

Wobei hier darauf zu achten ist dass 180°=π auf der Quincunx-Karte dem elliptischen Integral erster Art des Arguments 1/2 entspricht.

Sinusoidal

Code:

Eq = Import["C:\\Users\\Yukterez\\Desktop\\Equirectangular-projection.jpg"
(* Rohmaterial importieren *)]

Die Transformationen für die Sinusoidal-Projektion lauten hin und zurück:

Code:

Eq2Sin[{x_, y_}] := {(x - 0) Cos[y], y}
(* Sinusoidial <> Equirectangular *)
Sin2Eq[{x_, y_}] := {0 + x/Cos[y], y}
(* Equirectangular <> Sinusoidial *)

Wenn ein anderer Längengrad zentriert werden soll kann dieser in Radianten an der Stelle wo die 0 steht einfügt werden.
Transformation von der Plattkarte zum Sinusoidial:

Code:

Sn = ImageForwardTransformation[Eq, Eq2Sin,
  DataRange -> {{-π, π}, {-π/2, π/2}}, PlotRange -> {{-π, π}, {-π/2, π/2}}]

Code:

Sn = ImageTransformation[Eq, Sin2Eq,
  DataRange -> {{-π, π}, {-π/2, π/2}}, PlotRange -> {{-π, π}, {-π/2, π/2}}]
  (* Equirectangular -> Sinusoidial *)

Und vom Sinusoidial zurück in die Plattkarte:

Code:

ες = ImageTransformation[Sn, Eq2Sin,
  DataRange -> {{-π, π}, {-π/2, π/2}}, PlotRange -> {{-π, π}, {-π/2, π/2}}]

Code:

ες = ImageForwardTransformation[Sn, Sin2Eq,
  DataRange -> {{-π, π}, {-π/2, π/2}}, PlotRange -> {{-π, π}, {-π/2, π/2}}]
  (* Sinusoidial -> Equirectangular *)

Polar

Wie üblich zuerst die Transformation der Plattkarte in die mittabstandstreue Azimutalprojektion:

Code:

Eq = Import["C:\\Users\\Yukterez\\Desktop\\Equirectangular-projection.jpg"]
(* Plattkarte laden *)

Nachdem das Rohmaterial geladen ist die Funktion definieren:

Code:

Eq2Pl[{x_, y_}] := {(-2 y Cos[ x]), (2 y Sin[x])};
(* Polar Süd <> Equirectangular *)

Danach die als Variable Eq geladene Plattkarte in die Polarprojektion Pl transformieren:

Code:

Pl = ImageForwardTransformation[Eq, Eq2Pl,
  DataRange -> {{-π, π}, {0, π}},
  PlotRange -> {{-2 π, 2 π}, {-2 π, 2 π}}]
  (* Equirectangular -> Polar Süd *)

Falls sich bei zu niedrig aufgelöstem Rohmaterial an der Stelle wo sich die Ränder treffen eine Linie bildet kann man dieses Interpolationsartefakt mit einer leichten Modifikation der DataRange -> {{-π-Δ, π+Δ} in den Griff bekommen, wobei Δ=π/W und W für die Breite der Plattkarte in Pixeln steht (wenn das Rohmaterial 1000 px breit dann Δ=π/1000).

Liegt im umgekehrten Fall als Rohmaterial eine Polarprojektion die in die Rektangularprojektion transformiert werden soll vor wird die ImageForwardTransformation[] in eine ImageTransformation[] überführt:

Code:

εp = ImageTransformation[Pl, Eq2Pl,
  DataRange -> {{-2 π, 2 π}, {-2 π, 2 π}},
  PlotRange -> {{-π, π}, {0, π}}]
  (* Polar Süd -> Equirectangular *)

Soll statt des Südpols der Nordpol zentriert werden ändern wir das Vorzeichen:

Code:

Eq2P2[{x_, y_}] := {(2 y Cos[ x]), (2 y Sin[x])};
 (* Polar Nord <> Equirectangular *)

Und die DataRange:

Code:

pL = ImageForwardTransformation[Eq, Eq2P2,
  DataRange -> {{-π, π}, {-π, 0}},
  PlotRange -> {{-2 π, 2 π}, {-2 π, 2 π}}]
  (* Equirectangular -> Polar Nord *)

Daraus folgt im Umkehrschluss bei einer Rücktransformation

Code:

εP = ImageTransformation[pL, Eq2P2,
  DataRange -> {{-2 π, 2 π}, {-2 π, 2 π}},
  PlotRange -> {{-π, π}, {-π, 0}}]
  (* Polar Nord -> Equirectangular *)

Alternativ dazu kann man die Plattkarte vor der Transformation oder nach der Rücktransformation auch um 180° rotieren oder spiegeln wenn man für Nord- und Südzentrierung die selbe Transformation benutzen will oder irrtümlicherweise die falsche erwischt hat.

Stereographische Projektion

Zuerst wird eine Plattkarte in die konforme azimutale Projektion transformiert, und hernach wieder zurück. Mit dieser Projektion wird eine Kugelinversion, also die Invertierung einer Kugel in eine Kuppel, bewerkstelligt.

Code:

Eq = Import["C:\\Users\\Yukterez\\Desktop\\Equirectangular-projection.jpg"];
(* Plattkarte von der Festplatte laden *)

Die Funktionen lauten:

Code:

St2Eq[{x_, y_}] := {2 Tan[π/4 - y/2] Sin[x - 0], -2 Tan[π/4 - y/2] Cos[x - 0]};
Eq2St[{x_, y_}] := {0 + ArcTan[-y, x], -2 ArcTan[2, Sqrt[x^2 + y^2]] + π/2};
(* Equirectangular <> Stereographisch *)

Transformation Plattkarte Eq in die stereographische Projektion St:

Code:

Δ = 0;
St = ImageTransformation[Eq, Eq2St,
  DataRange -> {{-π - Δ, π + Δ}, {-π/2, π/2}}, PlotRange -> {{-2 π, 2 π}, {-2 π, 2 π}}, Padding-> "Fixed"]
(* Equirectangular -> Stereographisch *)

Stereographische Projektion St in eine Plattkarte eqSt zurücktransformieren:

Code:

Δ = 0;
eqSt = ImageTransformation[St, St2Eq,
  DataRange -> {{-2 π, 2 π}, {-2 π, 2 π}},  PlotRange -> {{-π - Δ, π + Δ}, {-π/2, π/2}}]
(* Stereographisch -> Equirectangular *)

Mercator

Nun wird eine Plattkarte in die Mercator-Projektion umgewandelt:

Code:

Eq = Import["C:\\Users\\Yukterez\\Desktop\\Equirectangular-projection.jpg"];
(* Plattkarte von der Festplatte laden *)

Die Transformationen:

Code:

Eq2Mc[{x_, y_}] := {x - 0, Log[Tan[π/4 + y/2]]};
Mc2Eq[{x_, y_}] := {x + 0, π/2 - 2 ArcTan[E^-y]};
(* Equirectangular <> Mercator *)

Transformation Plattkarte Eq in die Mercator Projektion Mc:

Code:

Mc = ImageForwardTransformation[Eq, Eq2Mc,
  DataRange -> {{-π, π}, {-π/2, π/2}},  PlotRange -> {{- π, π}, {- π, π}}]
(* Equirectangular \[Rule] Mercator *)

Mercator-Projektion Eq in eine Plattkarte eqMc zurücktransformieren:

Code:

eqMc = ImageForwardTransformation[Mc, Mc2Eq,
  DataRange -> {{-π, π}, {-π, π}}, PlotRange -> {{- π, π}, {- π/2, π/2}}]
(* Mercator -> Equirectangular *)

Karte neu zentrieren: Winkeltransformationen

Auf einer Plattkarte entspricht der Längengrad λ normalerweise der x-Achse, und der Breitengrad θ der y-Achse. So ergibt sich das bekannte Bild mit den Polen an den oberen und unteren Rändern und dem Äquator in der Mitte. Soll stattdessen eine andere Koordinate als {0,0} zentriert werden, z.B. der Nordpol, was einer Rotation um β=π/2=90° auf der θ-Achse entspricht, müssen die Kugelkoordinaten zuerst in den kartesischen 3D-Raum überführt, rotiert, und zurück in Kugelkoordinaten transformiert werden. Dazu schreiben wir zuerst folgende Funktion, um die Perspektive um die beliebigen Winkel α, β & ψ zu verschieben:

Code:

xyz[{x_, y_}] :=
 {
  Sin[y] Cos[x],
  Sin[y] Sin[x],
  Cos[y]
  }

Xyz[{x_, y_, z_}, α_] :=
 {
  x Cos[α] - y Sin[α],
  x Sin[α] + y Cos[α],
  z
  }

xYz[{x_, y_, z_}, β_] :=
 {
  x Cos[β] + z Sin[β],
  y,
  z Cos[β] - x Sin[β]
  }

xyZ[{x_, y_, z_}, ψ_] :=
 {
  x,
  y Cos[ψ] - z Sin[ψ],
  y Sin[ ψ] + z Cos[ψ]
  }

xy[{x_, y_, z_}] :=
 {
  ArcTan[x, y],
  ArcCos[z]
  }

rm[pic_, α_, β_, ψ_] :=
 xy[xyZ[xYz[Xyz[xyz[pic], α], β], ψ]]

RM[{x_, y_}] := rm[{x, y}, α, β, ψ]

Dann wird das Rohmaterial geladen:

Code:

Eq = Import["C:\\Users\\Yukterez\\Desktop\\Equirectangular-projection.jpg"]

Hier wird die Perspektive auf der x-Achse um den Winkel α gedreht:

Code:

α = +5/6 π; β = 0; ψ = 0;
p1 = ImageTransformation[Eq, RM,
  DataRange -> {{-π, π}, {0, π}},  PlotRange -> {{-π, π}, {0, π}}]

Jetzt auf der y-Achse um den Winkel β=90°, so dass der Nordpol zentriert ist:

Code:

α = 0; β = +π/2; ψ = 0;
p2 = ImageTransformation[Eq, RM,
  DataRange -> {{-π, π}, {0, π}},  PlotRange -> {{-π, π}, {0, π}}]

Nun mit β=-90°, also zentriertem Südpol:

Code:

α = 0; β = -π/2; ψ = 0;
p2 = ImageTransformation[Eq, RM,
  DataRange -> {{-π, π}, {0, π}},  PlotRange -> {{-π, π}, {0, π}}]

Im nächsten Bild wird der Winkel ψ (entlang der z-Achse) um 90° gedreht, also Norden und Süden nach Westen und Osten verschoben:

Code:

α = 0; β = 0; ψ = +π/2;
p3 = ImageTransformation[Eq, RM,
  DataRange -> {{-π, π}, {0, π}},  PlotRange -> {{-π, π}, {0, π}}]

Und eine ganze Umdrehung auf der y-Achse (β = -π..π):

Code:

α = 0; ψ = 0;
Do[Print[
  Eq2Eq[{x_, y_}] := rm[{x, y}, α, β, ψ];
  ImageTransformation[Eq, Eq2Eq,
  DataRange -> {{-π, π}, {0, π}}, PlotRange -> {{-π, π}, {0, π}}]],
 {β, -π, π, π/20}]

Wenn die selbe Operation an einer Mollweide, Hammer, Quincunx, Polar oder sonstigen Karte durchgeführt werden soll ist es am einfachsten diese wie in den vorangegangenen Beiträgen gezeigt zuerst in eine Plattkarte zu verwandeln, dann die neue Perspektive zu wählen und zuletzt wieder in die gewünschte Projektion zurückzutransformieren.

Sphärische 3D Projektion

Um eine Karte auf eine Kugel zu projizieren wandelt man diese wenn sie nicht bereits eine ist in eine Plattkarte um (siehe oben), da in dem Format die x- und y-Werte des Bildes 1:1 den Längen- und Breitengraden auf der Kugel entsprechen. Dann wird dieses Bild, hier Eq, als Textur definiert und auf die Kugel projiziert:

Code:

(* Methode 1 *)

image = Import["https://upload.wikimedia.org/wikipedia/commons/e/ea/Equirectangular-projection.jpg"];
Clear[α, β, ψ];
sphere = SphericalPlot3D[
 1, {θ, 0, π}, {φ, 0, 2 π},
 Mesh -> None, TextureCoordinateFunction -> ({#5, 1 - #4} &),
 PlotStyle -> Directive[Texture[image]],
 SphericalRegion -> True,
 Lighting -> "Neutral",
 Axes -> False,
 Boxed -> False,
 ViewPoint -> {100, 0, 20},
 ViewAngle -> 2/5,
 ImageSize -> 600,
 PlotPoints -> 150]
 
 
(* Methode 2 *)
 
image = Import["https://upload.wikimedia.org/wikipedia/commons/e/ea/Equirectangular-projection.jpg"];
Clear[α, β, ψ];
Kugel[pic_, {X_, Y_, Z_}] :=
 SphericalPlot3D[
 1, {u, 0, π}, {v, 0, 2 π},
 Mesh -> None,
 TextureCoordinateFunction -> ({#5, 1 - #4} &),
 PerformanceGoal -> "Quality",
 PlotStyle -> Directive[Texture[pic]],
 Lighting -> "Neutral",
 Axes -> False,
 RotationAction -> "Clip",
 SphericalRegion -> True,
 Boxed -> False,
 ViewPoint -> {X, Y, Z},
 ViewAngle -> 2/5,
 ImageSize -> 600,
 PlotPoints -> 150];
 sphere = Manipulate[Kugel[image, {
 +x Cos[α] Cos[β],
 +x Cos[ψ] Sin[α] + x Cos[α] Sin[β] Sin[ψ],
 -x Cos[α] Cos[ψ] Sin[β] + x Sin[α] Sin[ψ]}],
 {α, 0, 2 π}, {β, 0, 2 π}, {ψ, 0, 2 π}, {{x, 5}, 1/2, 20}]

Eine Übersicht der verschiedenen Projektionsarten, ihrer Eigenschaften und der sie beschreibenden Gleichungen findet sich hier und hier. Der gesamte Code für alle in diesem Faden verwendeten Transformationen in einer Datei befindet sich als .nb-File auf kartographie.yukterez.net

Für den Fall dass ein Gitter benötigt wird kann dieses wie folgt erzeugt:

Code:

color = Black; dX = H/9; dY = dX; W = 4000; H = W/2;
linesX := {color, Thick, Line[{{0 + dx, 0}, {0 + dx, H}}]};
linesY := {color, Thick, Line[{{0, 0 + dy}, {W, 0 + dy}}]};
lines = {Table[linesX, {dx, 0, W, dX}], Table[linesY, {dy, 0, H, dY}]};
grid = Graphics[lines, PlotRange -> {{0, W}, {0, H}}, Frame -> False]

und über das zu transformierende Bild Eq gelegt werden:

Code:

EQ = Show[Eq, grid]

H ist die Höhe des Rohmaterials Eq in Pixeln, und W dessen Breite. Alle Codes sind in der Wolfram-Sprache verfasst; zu den Anwendungsbeispielen und dem Grundgerüst in alternativer Syntax geht es hier entlang.
Cartographic Projections with Wolfram Mathematica. Formeln Gleichungen Equations