Kartographische Projektionen und Transformationen

[Yukterez]

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: Alles auswählen

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: Alles auswählen

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: Alles auswählen

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: Alles auswählen

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.image transformation bildtransformation bild mathematica geographie

[Yukterez]

Hammer-Aitoff

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

Code: Alles auswählen

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: Alles auswählen

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: Alles auswählen

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: Alles auswählen

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: Alles auswählen

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

hammer aitov hammeraitov hammeraitoff transformation equation formula projection

[Yukterez]

Quincunx

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

Code: Alles auswählen

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

Danach wird die Transformationsregel implementiert:

Code: Alles auswählen

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: Alles auswählen

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: Alles auswählen

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: Alles auswählen

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

Zurücktransformiert wird mit

Code: Alles auswählen

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. peirce quincuncial cartography transformation kartografie quintux

[Yukterez]

Sinusoidal

Code: Alles auswählen

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

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

Code: Alles auswählen

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: Alles auswählen

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

Code: Alles auswählen

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

Und vom Sinusoidial zurück in die Plattkarte:

Code: Alles auswählen

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

Code: Alles auswählen

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

Map Transformations with Mathematica

[Yukterez]

Polar

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

Code: Alles auswählen

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

Nachdem das Rohmaterial geladen ist die Funktion definieren:

Code: Alles auswählen

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: Alles auswählen

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: Alles auswählen

ε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: Alles auswählen

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

Und die DataRange:

Code: Alles auswählen

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

Daraus folgt im Umkehrschluss bei einer Rücktransformation

Code: Alles auswählen

ε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. Azimuthal equidistant projection transform equation world map

[Yukterez]

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: Alles auswählen

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

Die Funktionen lauten:

Code: Alles auswählen

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: Alles auswählen

Δ = 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: Alles auswählen

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

[Yukterez]

Mercator

Nun wird eine Plattkarte in die Mercator-Projektion umgewandelt:

Code: Alles auswählen

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

Die Transformationen:

Code: Alles auswählen

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: Alles auswählen

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

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

Code: Alles auswählen

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

[Yukterez]

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: Alles auswählen

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: Alles auswählen

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

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

Code: Alles auswählen

α = +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: Alles auswählen

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

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

Code: Alles auswählen

α = 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: Alles auswählen

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

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

Code: Alles auswählen

α = 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. transverse base projection, center equirectangular map focus angle latitude longitude 3d to 3d to 2d rotation Winkeltransformation

[Yukterez]

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: Alles auswählen

(* 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: Alles auswählen

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: Alles auswählen

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