Erzeugen von dreidimensionalen Graphiken mit scatterplot3d

Einordnung des Artikels

• Einführung in die Informatik
  • Einführung in die Programmiersprache R
    • Graphik
      • Erzeugen von dreidimensionalen Graphiken mit scatterplot3d

Da sich scatterplot3d() besonders einfach konfigurieren lässt, wenn die darzustellende Punktwolke als Dataframe gegeben ist, werden hier die Kenntnisse über Dataframes vorausgesetzt. Diese finden sich in Dataframes in R: der Datentyp data frame und Dataframes in R: Anwendungen.

Der Rückgabewert von scatterplot3d() ist eine Liste, deren Komponenten wiederum Funktionen sind; um diese aufzurufen, muss man mit dem Zugriff auf die Komponenten einer Liste vertraut sein. Alles Wissenswerte über Listen findet man in Listen in R: der Datentyp list und Listen in R: Anwendungen.

Weitere Beispiele, die das Paket scatterplot3d einsetzen (insbesondere die Funktion contour3d() zur Darstellung von Konturlinien), finden sich in Krummlinige Koordinatensysteme: Zylinderkoordinaten und Kugelkoordinaten. Anwendungen zum Zeichnen von Histogrammen mit xyz.convert() und segments() finden sich in Funktionsfabriken in R.

Einführung

Das Paket scatterplot3d ist eigentlich dazu gedacht Punktwolken darzustellen. Es besitzt aber derart reichhaltige Funktionalitäten, dass sich damit auch andere Aufgaben lösen lassen. Dabei sollte man aber bedenken, dass es dafür womöglich einfachere Möglichkeiten gibt.

Dieser Artikel soll nicht alle Funktionalitäten von scatterplot3d beschreiben, es werden lediglich einige besonders wichtige Aufgaben an Beispielen vorgestellt. Zumindest soll ihre ausführliche Besprechung dazu dienen, sich leicht in die Dokumentation einzuarbeiten.

Das zentrale Beispiel hier ist die Darstellung eines dreidimensionalen kubischen Gitters (siehe etwa Abbildung 1). Da die Gitterpunkte "in der Luft hängen", ist Abbildung 1 nicht geeignet, den richtigen dreidimensionalen Eindruck zu vermitteln; es werden einige Verbesserung vorgestellt, wie man dem Betrachter der Abbildung den Zugang zur dreidimensionalen Interpretation erleichtern kann und wie man weitere Informationen zu einem bestehenden Plot hinzufügen kann.

Das Paket scatterplot3d ist nicht in den Standardpaketen enthalten und muss erst installiert werden. Dies geschieht mit dem Befehl in Zeile 2 unten.

Um es in einem Skript einzusetzen, muss muss man den Befehl zum Laden des Paketes voranstellen (Zeile 4).

# Installation:
install.packages("scatterplot3d")
# Laden:
library("scatterplot3d")

Die erste Problemstellung: Darstellung eines kubischen Gitters

Es soll ein kubisches Gitter dargestellt werden, das aus einer Menge von Punkten (x, y, z) besteht, wobei

x = -5, -4, ..., 0, ..., 4, 5.

y = -5, -4, ..., 0, ..., 4, 5.

z = -1, 0, 1.

Das Gitter besteht somit aus 11·11·3 = 363 Gitterpunkten. Die Koordinaten der Gitterpunkte bilden die Testdaten, die später auf verschiedene Arten mit Hilfe von scatterplot3d dargestellt werden sollen.

Erzeugt werden die Testdaten – hier genannt coords –, indem man alle möglichen Kombinationen der angegebenen x-, y- und z-Werte bildet. Dies erledigt die Funktion expand.grid(). Diese Testdaten werden als ein Dataframe erzeugt (das wie eine Matrix eingesetzt werden kann):

n <- 5
xs <- (-n:n)
ys <- (-n:n)
zs <- (-1:1)
coords <- expand.grid(x = xs, y = ys, z = zs)

Das folgende Skript untersucht coords:

str(coords)
# 'data.frame': 363 obs. of  3 variables:
#   $ x: int  -5 -4 -3 -2 -1 0 1 2 3 4 ...
# $ y: int  -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 ...
# $ z: int  -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
# - attr(*, "out.attrs")=List of 2
# ..$ dim     : Named int  11 11 3
# .. ..- attr(*, "names")= chr  "x" "y" "z"
# ..$ dimnames:List of 3
# .. ..$ x: chr  "x=-5" "x=-4" "x=-3" "x=-2" ...
# .. ..$ y: chr  "y=-5" "y=-4" "y=-3" "y=-2" ...
# .. ..$ z: chr  "z=-1" "z= 0" "z= 1"

is.data.frame(coords)    # TRUE
is.matrix(coords)        # FALSE

str(coords$x)
# int [1:363] -5 -4 -3 -2 -1 0 1 2 3 4 ...

Das Dataframe coords besitzt somit 363 Zeilen und 3 Spalten; jede Zeile steht für die drei Koordinaten eines Gitterpunktes. Da beim Erzeugen des Dataframes (siehe erstes Skript Zeile 5) die Namen x, y und z für die Spalten gesetzt wurden, kann man die Spalten

• entweder wie bei einer Matrix, etwa mit coords[ , 1] für die erste Spalte,
• oder mit dem Namen, etwa coords$x , auswählen.

Darstellung des Gitters mit scatterplot3d

Die einfache Darstellung

Das folgende Skript stellt die 363 Gitterpunkte im dreidimensionalen Koordinatensystem dar, dazu wird die Funktion scatterplot3d() mit ihren wichtigsten Argumenten eingesetzt; das Ergebnis ist in Abbildung 1 zu sehen:

scatterplot3d(x = coords[ , 1], y = coords[ , 2], z = coords[ , 3], 
                    color = "red", type = "p",
                    xlab = "x", ylab = "y", zlab = "z",
                    lwd = 3, angle = 40,
                    main = "Gitter")

Zur Erklärung:

Zeile 1: Die x-, y- und z-Koordinaten werden mit den Argumenten x, y und z gesetzt. Da coords 363 Zeilen hat werden somit alle 363 Gitterpunkte gezeichnet. Suggestiver und vielleicht weniger fehleranfällig ist die Schreibweise scatterplot3d(x = coords$x, y = coords$y, z = coords$z) zum Setzen der Koordinaten. Für das Ergebnis sind beide Schreibweisen gleichwertig, Letztere beruht darauf, dass die Namen der Spalten gesetzt sind. Dass man anstelle der drei Koordinaten x, y und z viel einfacher das Dataframe coords übergeben kann, wird weiter unten gezeigt.

Zeile 2: Mit color = "red" wird die Farbe der Gitterpunkte gesetzt, mit type = "p&" werden tatsächlich Punkte erzeugt. Man beachte, dass die Farbe hier nicht wie bei der Funktion plot() mit col = "red" gesetzt wird. Und type = "p" ist default-Wert, den man nicht setzen muss.

Zeile 3: Die Achsen werden beschriftet, indem man xlab = "x", ylab = "y", zlab = "z" setzt.

Zeile 4: Die Strichdicke wird mit lwd = 3 gesetzt. Der Winkel zwischen der x- und y-Achse, der für die perspektivische Darstellung sorgt, wird durch angle = 40 definiert. Die Angabe von 40 (Grad im Winkelmaß) wäre nicht nötig gewesen, da dies der default-Wert ist. Je nachdem, wie die darzustellenden Objekte liegen, wird man den Winkel geeignet anpassen.

Zeile 5: Mit main = "Gitter" wird der Titel der Abbildung gesetzt.

Im Skript wurden mehrere Größen nicht gesetzt; dass sie dennoch gezeigt sind, beruht darauf, dass entsprechende default-Werte verwendet werden:

• Der Wertebereich für x, y und z, in dem das Koordinatensystem gezeigt wird, wird automatisch aus den darzustellenden Werten berechnet. Man könnte ihn mit xlim, ylim und zlim anpassen.
• Dass die Achsen und ihre Beschriftung gezeichnet werden, beruht auf den default-Werten axis=TRUE, tick.marks=TRUE, label.tick.marks=TRUE .
• Dass eine Box um den Plot und ein Gitternetz auf dessen Boden gezeichnet wird, beruht auf den default-Werten grid=TRUE, box=TRUE

Abbildung 1: Darstellung des kubischen Gitters durch Punkte. Viele Graphik-Parameter der Funktion scatterplot3d() sind durch default-Werte sinnvoll gesetzt und müssen daher nicht selbst konfiguriert werden.

Darstellung als Histogramm

Bei vielen Anwendungen, in denen Punktwolken vorkommen, ist ein Histogramm besser geeignet als die Darstellung von Punkten. Der Nachteil von Punkten ist häufig, dass sie "in der Luft hängen" und daher die dreidimensionale Interpretation sehr schwierig oder nicht eindeutig ist. Manchmal lassen sich Histogramme besser lesen. Dazu wird lediglich in obigem Skript type = "p" durch type = "h" ersetzt. Zusätzlich wurde lwd = 2 gesetzt. Das Ergebnis zeigt Abbildung 2; die Veränderungen führen hier aber nicht zu einer leichter erkennbaren Darstellung des Gitters – bei anderen Punktmengen mögen Histogramme besser geeignet sein.

scatterplot3d(x = coords[ , 1], y = coords[ , 2], z = coords[ , 3], 
                    color = "red", type = "h",
                    xlab = "x", ylab = "y", zlab = "z",
                    lwd = 2, angle = 40,
                    main = "Gitter")

Abbildung 2: Im Vergleich zur Darstellung des kubischen Gitters in Abbildung 1 sind zwei Graphik-Parameter in scatterplot3d() verändert: Anstelle der Punkte werden Histogramme gezeichnet (was in diesem speziellen Beispiel keine Verbesserung darstellt) und die Linienbreite lwd (linewidth) wurde von drei auf zwei Pixel verringert.

Die Übergabe der Koordinaten an scatterplot3d()

Die zentralen Eingabewerte von scatterplot3d() sind natürlich x, y und z, mit denen die Koordinaten der Punkte in der Wolke gesetzt werden. In der Dokumentation wird angegeben (weitere Argumente sind hier weggelassen):

scatterplot3d(x, y=NULL, z=NULL)

Auf den ersten Blick mag hier verwundern, warum y und z per default mit NULL gesetzt werden. Dies soll andeuten, dass sie nicht gesetzt werden müssen, sondern dass man die Koordinaten in einem geeigneten Objekt x verpacken kann. Und dieses Objekt ist natürlich ein Dataframe, in dem die darzustellenden Daten oft abgelegt sind. (Umgekehrt ist es meist kein großer Aufwand, die Daten entsprechend aufzubereiten.)

Anstelle des ersten Skripts kann man daher einfacher schreiben:

scatterplot3d(x = coords, 
                    color = "red", type = "p",
                    xlab = "x", ylab = "y", zlab = "z",
                    lwd = 3, angle = 40,
                    main = "Gitter")

Im folgenden wird stets diese Schreibweise verwendet und coords verwendet, das im ersten Skript mit expand.grid() erzeugt wurde.

Verbesserte dreidimensionale Darstellung

Ein Problem für die dreidimensionale Interpretation eines Plots wie Abbildung 1 ist, dass sich die Punkte überlagern. Hier sind etwa Punkte aus Ebenen zu verschiedenen z-Werten nicht klar voneinander zu trennen. Oft – es gibt natürlich auch Gegenbeispiele – reicht es dazu, den Winkel angle zwischen x- und y-Achse anzupassen. Das folgende Skript verwendet angle = 27 anstelle des default-Wertes von 40 (siehe Zeile 4); das Ergebnis ist in Abbildung 3 gezeigt:

scatterplot3d(x = coords, 
                color = "red", type = "p",
                xlab = "x", ylab = "y", zlab = "z",
                lwd = 3, angle = 27,
                main = "Gitter")

Abbildung 3: Im Vergleich zur Darstellung des kubischen Gitters in Abbildung 1 ist der Winkel angle in scatterplot3d() verändert: Durch die Wahl von 27 Grad wird vermieden, dass sich die Punkte von Ebenen mit unterschiedlichen z-Werten überlagern.

Eine weitere Verbesserung lässt sich manchmal erzielen, indem man highlight.3d = TRUE setzt (der default-Wert ist FALSE; siehe Zeile 5). Jetzt werden die Punkte entlang der y-Achse unterschiedlich gefärbt. Das Ergebnis zeigt Abbildung 4.

scatterplot3d(x = coords, 
                    color = "red", type = "p",
                    xlab = "x", ylab = "y", zlab = "z",
                    lwd = 3, angle = 27,
                    highlight.3d = TRUE,
                    main = "Gitter")

Abbildung 4: Durch Setzten von highlight.3d = TRUE in scatterplot3d() werden die Punkte gemäß ihrer y-Koordinate unterschiedlich gefärbt, wodurch der räumliche Eindruck verstärkt wird.

Gruppierung der Punkte und ihre Kennzeichnung durch Farben

In vielen Anwendungen ist eine Gruppierung der Punktmenge inhaltlich sinnvoll. Diese Gruppierung kann durch unterschiedliche Farben sichtbar gemacht werden.

Im vorliegenden Beispiel sind die x- und y-Koordinaten austauschbar (beide laufen von -5 bis 5), dagegen sind die z-Koordinaten ausgezeichnet, da sie lediglich drei Werte annehmen können. Daher ist eine Gruppierung (oder Schichtung) gemäß der z-Koordinate sinnvoll.

Dazu werden drei unterschiedliche Farben vereinbart, die hier mit Hilfe der Palette rainbow gebildet werden. (Mehr Informationen über Farbpaletten findet man im Paket grDevices unter Palettes.) Im folgenden Skript wird das Dataframe coords ausgewertet und bestimmt, wie viele unterschiedliche z-Koordinaten (Schichten) enthalten sind (es wird also nicht verwendet, dass es drei unterschiedliche z-Werte gibt, siehe Zeile 2 und 4). Entsprechend wird die Anzahl der Farben in der Farbpalette gewählt (Zeile 6). Anschließend wird ein Vektor aus Farben gebildet, dessen Länge mit der Anzahl der Zeilen von coords übereinstimmt (Zeile 10); die unterschiedlichen z-Werte werden jeweils einer Farbe zugeordnet (Zeile 11 bis 13). Die entsprechenden Anweisungen sind durch Kommentare erläutert.

# Anzahl der Zeilen von coords:
lgth <- nrow(coords)
# Anzahl der Schichten:
n.colors <- length(unique(coords[ , 3]))
# entspr. Anzahl der Farben:
palette.colors <- rainbow(n = n.colors)

# Färben gemäß z-Koord.:
z.min <- min(coords[ , 3])
colors <- vector(mode = "character", length = lgth)
for(i in (1:lgth)){
  colors[i] <- palette.colors[coords[i, 3] + 1 - z.min]
}

str(colors)
# chr [1:363] "#FF0000FF" "#FF0000FF" "#FF0000FF" "#FF0000FF" "#FF0000FF" "#FF0000FF" "#FF0000FF" "#FF0000FF" ...

Das Skript wirkt auf den ersten Blick sehr umständlich, denn es lässt sich aus den möglichen z-Werten zs <- (-1:1) (vor der Definition von coords) einfach ablesen, dass es drei unterschiedliche z-Werte gibt. Der Vorteil des Skripts liegt darin, dass durch length(unique(coords[ , 3])) die Anzahl der unterschiedlichen z-Werte in coords bestimmt wird (Zeile 4). Das heißt bei einer anderen Definition von coords kann dieser Befehl unverändert bleiben und er liefert wieder die Anzahl der unterschiedlichen z-Werte. Auch wie in Zeile 9 und 12 die Zuordnung zwischen den möglichen z-Werten und den Farben hergestellt wird, muss nicht angepasst werden. (Sie beruht allerdings darauf, dass sich die aufeinanderfolgenden z-Werte um jeweils 1 unterscheiden.)

Im Aufruf von scatterplot3d() wird jetzt nicht eine einzige Farbe gesetzt wie in der obigen Skripten, sondern es wird der Vektor colors verwendet. Da dessen Länge mit der Anzahl von Zeilen von coords übereinstimmt, erfolgt eine Zuordnung von jedem Gitterpunkt aus dem Dataframe coords zu der entsprechenden Farbe in colors; das Ergebnis ist in Abbildung 5 zu sehen.

scatterplot3d(x = coords, 
                    color = colors, type = "p",
                    xlab = "x", ylab = "y", zlab = "z",
                    lwd = 3, angle = 27,
                    main = "Gitter")

Abbildung 5: Im Aufruf von scatterplot3d() wird die Farbe color mit Hilfe des Vektors colors gesetzt, dessen Länge mit der Anzahl der Zeilen von coords übereinstimmt. Dadurch kann eigentlich jedem Gitterpunkt eine eigene Farbe zugeordnet werden; so wie hier colors definiert ist, werden aber die Farben den unterschiedlichen z-Koordinaten zugeordnet.

Im Skript oben wurden mehrere Variablen definiert, die eigentlich nicht benötigt werden, da man die Funktionsaufrufe geeignet verschachteln kann; es ist nur geschehen, um zu verdeutlichen, welche Berechnungen ausgeführt werden müssen, um die Farben zu setzen.

Das folgende Skript setzt für jeden Gitterpunkt eine eigene Farbe und vermeidet unnötige Variablen (lediglich das Dataframe coords wird als gegeben vorausgesetzt). Das Ergebnis ist in Abbildung 6 zu sehen.

scatterplot3d(x = coords, 
              color = rainbow(n = nrow(coords)), type = "p",
              xlab = "x", ylab = "y", zlab = "z",
              lwd = 3, angle = 27,
              main = "Gitter")

Abbildung 6: Darstellung des kubischen Gitters mit scatterplot3d(), so dass jeder Gitterpunkt eine eigene Farbe erhält; verwendet wird dazu wie in Abbildung 5 die Farbpalette rainbow.

Der Rückgabewert von scatterplot3d()

Der Rückgabewert als Liste von sechs Funktionen

In allen bisherigen Anwendungen wurde die Funktion scatterplot3d() in dem Sinn eingesetzt, dass sie Anweisungen ausführt, aber nicht in dem Sinn, dass sie einen Rückgabewert berechnet, der abgespeichert und weiter verwendet werden kann. Um zu testen, welchen Rückgabewert die Funktion scatterplot3d() erzeugt, wird ein Funktionsaufruf von oben ausgeführt, als Objekt sc abgespeichert und dessen Struktur untersucht:

sc <- scatterplot3d(x = coords, 
              color = "red", type = "p",
              xlab = "x", ylab = "y", zlab = "z",
              lwd = 3, angle = 27,
              main = "Gitter")

str(sc)
# List of 6
# $ xyz.convert:function (x, y = NULL, z = NULL)  
#   $ points3d   :function (x, y = NULL, z = NULL, type = "p", ...)  
#     $ plane3d    :function (Intercept, x.coef = NULL, y.coef = NULL, lty = "dashed", lty.box = NULL, draw_lines = TRUE, draw_polygon = FALSE, 
#                             polygon_args = list(border = NA, col = rgb(0, 0, 0, 0.2)), ...)  
#       $ box3d      :function (...)  
#         $ contour3d  :function (f, x.count = 10, y.count = 10, type = "l", lty = "24", x.resolution = 50, y.resolution = 50, ...)  
#           $ par.mar    :List of 1
# ..$ mar: num [1:4] 5.1 3.1 4.1 3.1

Man erkennt, dass eine Liste mit 6 Komponenten zurückgegeben wird; die Komponenten sind wiederum Funktionen, deren Eingabewerte man an der Ausgabe der Struktur im Skript oben ablesen kann. Die Aufgabe der 6 Funktionen wird in folgender Tabelle kurz beschrieben:

Es sollen hier nicht sämtliche Eigenschaften dieser 6 Funktionen vorgestellt werden; die Erklärungen beschränken sich auf einige Spezialfälle:

1. Hinzufügen von Ebenen zu einem gegebenen Plot mit plane3d().
2. Hinzufügen von Punkten mit points3d().
3. Hinzufügen von Linien zu einem Plot. Dazu wird hier die Funktion segments() aus dem Paket graphics verwendet. Um segments() einzusetzen, müssen aber die Projektionen der Endpunkte der Linien in die zweidimensionale Ebene berechnet werden, wozu die Funktion xyz.convert() verwendet wird.

Zum Zeichnen von Konturlinien mit contour3d() wird hier kein Beispiel gezeigt; in Krummlinige Koordinatensysteme: Zylinderkoordinaten und Kugelkoordinaten sind einige Anwendungen mit den Quelltexten zu finden.

Hinzufügen von Ebenen zu einem gegebenen Plot mit plane3d()

Dadurch dass die Gitterpunkte wie in Abbildung 5 und 6 gefärbt wurden, ist das Bild zwar leichter räumlich zu interpretieren, dennoch "schweben" die Gitterpunkte in der Luft, so dass die räumliche Interpretation nicht eindeutig ist. Man kann diesen Effekt ein wenig abschwächen, indem man die entsprechenden xy-Ebenen in den Plot einfügt – ähnlich wie auf dem Boden der umgebenden Box das Koordinatennetz gezeichnet ist.

Dass Ebenen in einen Plot eingetragen werden, geschieht natürlich nicht nur aus dem Grund der leichteren räumlichen Interpretation, oft haben diese Ebenen eine bestimmte inhaltliche Bedeutung und sind daher ein wichtiges Instrument, um dreidimensionale Abbildungen gehaltvoll zu gestalten.

Um zu verstehen, wie die Funktion plane3d() zum Zeichnen von Ebenen eingesetzt wird, muss man vor allem die drei Eingabewerte Intercept, x.coef = NULL, y.coef = NULL kennen, mit denen die Ebene beschrieben wird. Dazu wird die Ebene in der Form

z(x, y) = z₀ + a·x + b·y

dargestellt. Sofern die Ebene nicht parallel zur z-Achse ist, ist diese Darstellung möglich und die drei Parameter z₀, a und b sind eindeutig bestimmt; sie liefern drei Eingabewerte Intercept, x.coef und y.coef.

Das folgende Skript versucht Abbildung 5 zu verbessern, indem mit Hilfe plane3d() die drei Ebenen eingetragen werden, auf denen sich Gitterpunkte mit identischen z-Koordinaten befinden. Die Zeilen 1 bis 13 zeigen nochmals wie der Vektor colors der Farben der Gitterpunkte definiert wird (identisch wie zu Abbildung 5).

# Anzahl der Zeilen von coords:
lgth <- nrow(coords)
# Anzahl der Schichten:
n.colors <- length(unique(coords[ , 3]))
# entspr. Anzahl der Farben:
palette.colors <- rainbow(n = n.colors)

# Färben gemäß z-Koord.:
z.min <- min(coords[ , 3])
colors <- vector(mode = "character", length = lgth)
for(i in (1:lgth)){
  colors[i] <- palette.colors[coords[i, 3] + 1 - z.min]
}

scp <- scatterplot3d(x = coords, 
              color = colors, type = "p",
              xlab = "x", ylab = "y", zlab = "z",
              lwd = 3, angle = 27,
              main = "Gitter")

lapply( X = (1:n.colors), 
        FUN = function(i){scp$plane3d(Intercept = z.min + i - 1, x.coef = 0, y.coef = 0,
                                      lty = 1,
                                      col = palette.colors[i]) } )

Zeile 15 bis 19: Der Aufruf von scatterplot3d() beinhaltet nur einen Unterschied im Vergleich zum Skript, mit dem Abbildung 5 erzeugt wurde – aber dieser Unterschied erlaubt den Aufruf der oben kurz vorgestellten 6 Funktionen: Der Rückgabewert wird in der Variable scp abgespeichert (scp soll an scatterplot3d erinnern, Zeile 15); dies ist eine Liste von Funktionen, die später verwendet wird, um die Ebenen hinzuzufügen. Ansonsten sind die Argumente in scatterplot3d() wie in den besprochenen Skripten gesetzt.

Zeile 21 bis 24: Wirklich neu ist jetzt der Aufruf von plane3d(), der innerhalb der Schleife geschieht, die mit lapply() realisiert wird (Zeile 21). Dabei wird über die unterschiedlichen z-Werte der Ebenen von Gitterpunkten (beziehungsweise deren Farben) iteriert, hier also über drei Ebenen oder Farben. In jedem Iterationsschritt wird das Objekt scp verwendet, um die Funktion plane3d() aufzurufen (Zeile 22). Innerhalb plane3d() wird die Ebene durch die drei Argumente Intercept, x.coef und y.coef definiert. An Zeile 9 und 12 erkennt man, wie die Zuordnung zwischen der z-Koordinate und der Farbe hergestellt wird (der kleinste z-Wert gehört zur ersten Farbe der Palette, also rot bei rainbow, der größte z-Wert zur letzten Farbe der Palette, also blau). Da alle Ebenen parallel zur xy-Ebene sind, sind x.coef = 0 und y.coef = 0 . Mit lty = 1 (lty = linetype) wird festgelegt, dass die Linien durchgezogen und nicht gestrichelt (default-Wert) gezeichnet werden. Und mit col = palette.colors[i] wird die Farbe der Ebene festgelegt. Man beachte, dass hier das Argument col heißt und nicht color.

Abbildung 7 zeigt das Ergebnis des Skripts, das man jetzt mit Abbildung 1 vergleichen kann, um zu beurteilen, ob die zusätzlichen Maßnahmen zu einer Verbesserung der Darstellung geführt haben.

Abbildung 7: Darstellung des kubischen Gitters mit scatterplot3d(), wobei die Gitterpunkte auf jeder der drei Ebenen mit identischen z-Werten eine eigene Farbe erhalten; verwendet wird dazu wie in Abbildung 5 die Farbpalette rainbow. Zusätzlich werden mit Hilfe der Funktion plane3d() die drei Ebenen eingetragen, wodurch der Eindruck verstärkt wird, dass sich die gleichfarbigen Punkte tatsächlich auf einer gemeinsamen Ebene befinden.

Hinzufügen von Punkten zu einem gegebenen Plot mit points3d() und die Funktion xyz.convert()

Es sollen nun der Plot aus Abbildung 7 ergänzt werden:

1. Es werden drei Punkte hinzugefügt und als Histogramm eingetragen (schwarz und gestrichelt). Ihre Koordinaten sollen zufällig gewählt werden, sich aber innerhalb des Bereichs befinden, in dem auch das Gitter liegt.
2. Die drei Punkte werden durch Strecken miteinander verbunden (schwarz mit durchgezogenen Linien).

Das folgende Skript zerlegt diese Aufgaben in 4 Teilaufgaben, dabei sind die Vorbereitungen, also coords und den Vektor der Farben erzeugen, nicht nochmals aufgeführt:

1. Mit runif werden die Koordinaten der drei neuen Punkte berechnet; das Ergebnis wird in einem Dataframe randomPoints abgelegt.
2. Der Plot aus Abbildung 7 mit den drei Ebenen wird erzeugt. Der Rückgabewert als scp wird abgespeichert, um die weiteren Aufgaben mit den Funktionsaufrufen von points3d() und xyz.convert() zu erledigen.
3. Die zufällig erzeugten Punkte werden als Histogramm eingetragen. Hier wird die Funktion points3d() benötigt.
4. Die Verbindungslinien zwischen randomPoints eintragen: Dies wird sich als die schwierigste Aufgabe herausstellen, da zum Zeichnen von Linien eine Funktion wie segments() benötigt wird, die eigentlich für zweidimensionale Graphiken implementiert ist. Um sie mit geeigneten Koordinaten zu füttern, müssen die dreidimensionalen Koordinaten aus randomPoints erst mit Hilfe von xyz.convert() in die zweidimensionale Ebene projiziert werden.

# 1. Erzeugen der drei Punkte, Abspeichern im Dataframe randomPoints:
set.seed(seed = 7)
N <- 3
v.5.x <- runif(n = 3, min = -5, max = 5)
v.5.y <- runif(n = 3, min = -5, max = 5)
v.1.z <- runif(n = 3, min = -1, max = 1)
randomPoints <- as.data.frame( cbind(x = v.5.x, y = v.5.y, z = v.1.z) )
print(randomPoints, digits = 3)
#      x     y      z
# 1  4.89 -4.30 -0.320
# 2 -1.02 -2.56  0.944
# 3 -3.84  2.92 -0.668

# 2. Plot aus Abbildung 7 erzeugen:
scp <- scatterplot3d(x = coords,
                     color = colors, type = "p",
                     xlab = "x", ylab = "y", zlab = "z",
                     lwd = 3, angle = 27,
                     main = "Gitter")

# Gitterebenen eintragen:
lapply( X = (1:n.colors), 
        FUN = function(i){scp$plane3d(Intercept = z.min + i - 1, x.coef = 0, y.coef = 0,
                                      lty = 1,
                                      col = palette.colors[i]) } )

# 3. randomPoints als Histogramm eintragen:
scp$points3d(x = randomPoints,
             type = "h", lty = "dotted",
             lwd = 3, col = "black")

# 4. Projektion der randomPoints berechnen und 
# Verbindungslinien zwischen den Punkten eintragen: 
p.2 <- scp$xyz.convert(x = randomPoints)
str(p.2)
# List of 2
# $ x: num [1:3] 4.51 1.58 2.25
# $ y: num [1:3] 1.606 0.661 -0.564

lapply( X = (1:N), FUN = function(i)
  {
    if(i < N){j <- i+1} else {j <- 1};
    segments(x0 = p.2$x[i], y0 = p.2$y[i], x1 = p.2$x[j], y1 = p.2$y[j], lwd = 3)
  } )

Es werden nur die neuen Bestandteile des Skripts kurz erklärt.

Zeile 2: Damit das Skript reproduziert werden kann, wird der Zufallsgenerator explizit initialisiert.

Zeile 7: Die Koordinaten der randomPoints werden wieder in einem Dataframe abgelegt, da man dies in den Funktionen einfacher aufrufen kann als die einzelnen Koordinaten.

Zeile 15 bis 25: wurden bereits ausführlich besprochen.

Zeile 28 bis 30: Mit dem Objekt scp (siehe Zeile 15) wird die Funktion points3d() aufgerufen. Sie erhält als Eingabewerte das Dataframe randomPoints und weitere Graphik-Parameter.

Zeile 34: Ebenso benötigt man scp zum Aufruf von xyz.convert(). Sie erhält als Eingabewert das Dataframe randomPoints und berechnet daraus zweidimensionale Koordinaten der drei Punkte. Sie werden als p.2 abgespeichert (p soll an point erinnern, 2 an zweidimensionale Koordinaten).

Zeile 35 bis 38: An der Ausgabe der Struktur von p.2 erkennt man, wie xyz.convert() die zweidimensionalen Koordinaten zu einem Objekt zusammenfügt. Es entsteht eine Liste mit zwei Komponenten, je eine Komponente für die x- beziehungsweise y-Koordinaten. Und da drei Punkte in randomPoints enthalten waren, ist jede dieser Komponenten ein Vektor der Länge 3.

Zeile 43: Die Funktion segments() wird verwendet, um die Verbindungslinien der 3 Punkte zu zeichnen. Sie befindet sich im Paket graphics und besitzt 4 Koordinaten als Eingabewerte sowie weitere Graphik-Parameter: segments(x0, y0, x1 = x0, y1 = y0) . Die Graphik-Parameter sind hier nicht gezeigt, die 4 Koordinaten erklären sich selbst, wenn man die Aufgabe der Funktion kennt.

Zeile 40 bis 44: Innerhalb der Schleife, die hier mit lapply() realisiert ist, durchläuft man die drei Punkte und zeichnet die Verbindungslinien. In Zeile 43 erkennt man, wie p.2 eingesetzt wird, um die Argumente von segments() zu setzen.

Das Ergebnis des Skripts ist in Abbildung 8 dargestellt.

Abbildung 8: Zusätzlich zu dem Plot aus Abbildung 7 werden mit Hilfe der Funktionen points3d() und xyz.convert() Punkte und Verbindungslinien eingezeichnet. Die Koordinaten der Punkte wurden zuvor mit einem Zufallsgenerator erzeugt und in einem Dataframe abgespeichert. Die Funktion points3d() kann wieder leicht auf das Dataframe zugreifen und zeichnet die Punkte (wozu geeignete Graphik-Parameter gesetzt werden müssen). Die Funktion xyz.convert() berechnet die Projektion der Punkte in ein zweidimensionales Koordinatensystem, so dass die Funktion segments() verwendet werden kann, um die Verbindungslinien der Punkte zu zeichnen.

Aufgabe: Wie muss man das Skript abändern, um anstelle de schwarzen Dreiecks einen Polygonzug mit einer anderen Eckenzahl zu erzeugen. Führen Sie dies aus. (Achtung: bei großer Eckenzahl wird der Plot sehr unübersichtlich!)

Die zweite Problemstellung: Darstellung einer Irrfahrt (Random Walk)

Definition der Irrfahrt (Random Walk)

Es soll eine Irrfahrt (oder Random Walk) dargestellt, die durch folgende Bedingungen definiert ist:

1. Ein Teilchen kann sich auf einem kubischen (also dreidimensionalen) Gitter bewegen; das Gitter ist parallel zu den Koordinatenachsen und zwei benachbarte Gitterpunkte besitzen den Abstand 1.
2. Damit ist das Gitter aufgebaut wie das aus der ersten Problemstellung, mit dem Unterschied, dass es sich in alle Richtungen ins Unendliche erstreckt.
3. Das Teilchen startet im Ursprung (0, 0, 0) und kann sich in jedem Schritt um 1 in eine der drei Raumrichtungen bewegen.
4. Die Zufälligkeit der Bewegung ist durch zwei Elemente gegeben: Zum Einen wird die Richtung aus den drei Möglichkeiten (x-, y- oder z-Richtung) mit jeweils gleicher Wahrscheinlichkeit 1/3 ausgewählt. Zum Anderen bewegt es sich mit Wahrscheinlichkeit 0.55 jeweils in positive Richtung und nur mit Wahrscheinlichkeit 0.45 in negative Richtung. Auf lange Sicht wird sich das Teilchen mit hoher Wahrscheinlichkeit im ersten Oktanten aufhalten, also im Bereich mit x > 0, y > 0 und z > 0.
5. Betrachtet werden sollen 1000 Schritte des Teilchens. Der Quelltext soll aber so gestaltet werden, dass man diese Zahl leicht abändern kann.

Darstellung der Irrfahrt

Die folgende Abbildung 9 zeigt eine mögliche Realisierung der Irrfahrt. In der Graphik sind folgende Elemente enthalten:

1. Der Pfad des Teilchens ist rot eingetragen; man erkennt deutlich, dass in jedem Schritt eine Bewegung entlang einer der Koordinatenachsen stattfindet. Allerdings sind die Koordinatenachsen sehr unterschiedlich skaliert, was die Interpretation der Graphik erschwert.
2. Die xy-Ebene zu z = 0 ist schwarz und mit gestrichelten Linien eingetragen.
3. Um die räumliche Interpretation der Graphik zu erleichtern, sind die Punkte nach je einem Zeitschritt als Histogramm eingetragen; sie sind grün und gestrichelt.
4. Die Spur der Bewegung in der xy-Ebene, also die Projektion des Pfades auf die xy-Ebene (zu z = 0), ist blau eingetragen.
5. Zur besseren Orientierung sind der Ursprung als Startpunkt der Bewegung, die Gitterpunkte auf der Raumdiagonale, also die Punkte (0, 0, 0), (1, 1, 1) ..., sowie die Spur der Raumdiagonale in der xy-Ebene eingetragen (schwarz).

Abbildung 9: Die oben beschriebene Irrfahrt wird einmal mit 1000 Schritten realisiert und dargestellt. Zum Erzeugen der Graphik werden die Funktionen scatterplot3d(), plane3d(), points3d() sowie xyz.convert() verwendet. Welche Elemente außer dem Pfad des Teilchens dargestellt sind, wird im Text näher erklärt.

Das Skript zum Erzeugen der Irrfahrt und der Graphik

Unten wird das Skript gezeigt, mit dem eine Realisierung der Irrfahrt sowie ihre Darstellung in Abbildung 9 erzeugt wird. Das Skript kann in 9 Teilaufgaben unterteilt werden, die hier nur grob beschrieben werden:

1. Erzeugen des Dataframes path, das die Koordinaten des Teilchens während der 1000 Schritte enthält. Dabei wird die Irrfahrt so realisiert, wie sie oben beschriebenen wurde.
2. Darstellung des Pfades path mit scatterplot3d() (rot).
3. Eintragen von Ursprung und xy-Ebene in den Plot (schwarz).
4. Vorbereitungen für das Histogramm (grün) und die Spur in der xy-Ebene (blau).
5. Eintragen der Verbindungen zwischen dem Pfad path und seiner Spur in der xy-Ebene (grünes Histogramm).
6. Eintragen der Spur in der xy-Ebene, also Verbindung aufeinanderfolgender Punkte in der Spur von path (blau).
7. Raumdiagonale berechnen.
8. Raumdiagonale eintragen (schwarze Punkte).
9. Spur der Raumdiagonale eintragen (schwarz).

# 1. Realisierung der Irrfahrt: Erzeugen des Dataframes path
set.seed(seed = 42)
lgth <- 1000
path <- as.data.frame( matrix(data = 0, nrow = lgth, ncol = 3) )

# Richtung (positiv oder negativ):
deltas <- sample(x = c(-1, 1), size = lgth - 1, replace = TRUE, prob = c(0.45, 0.55))

# Richtung (x, y oder z):
idxs <- sample(x = (1:3), size = lgth - 1, replace = TRUE)

# Erzeugen von path (Start im Ursprung):
for( i in 1:(lgth-1) ){
  path[i+1, ] <- path[i, ]
  idx <- idxs[i]
  path[i+1, idx] <- path[i, idx] + deltas[i]
}

# str(path)
# 'data.frame': 1000 obs. of  3 variables:
#   $ V1: num  0 0 0 1 1 0 0 0 0 0 ...
# $ V2: num  0 0 -1 -1 -2 -2 -1 -2 -1 -2 ...
# $ V3: num  0 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...

# 2. Darstellung des Pfades path mit scatterplot3d():
scp <- scatterplot3d(x = path, 
                     color = "red",
                     type = "l",
                     xlab = "x", ylab = "y", zlab = "z",
                     lwd = 2, angle = 75,
                     main = "Random Walk")

# 3. Eintragen von Ursprung und xy-Ebene
# Ursprung:
scp$points3d(x = path[1, ],
             type = "p",
             lwd = 3, col = "black")
# xy-Ebene:
scp$plane3d( Intercept = 0, x.coef = 0, y.coef = 0, 
                 lty = "dotted", col = "black" )

# 4. Vorbereitungen für das Histogramm (grün) und die Spur in der xy-Ebene (blau):
# Projektion von path auf xy-Ebene (z-Komponente gleich 0 setzen)
path.xy <- cbind( path[ , c(1,2)], rep(x = 0, times = lgth) )

# Projektion nach 2D:
path.xy.2 <- scp$xyz.convert(x = path.xy)
path.2 <- scp$xyz.convert(x = path)

# 5. Eintragen der Verbindungen zwischen path und path.xy (grünes Histogramm)
lapply( X = (1:lgth), FUN = function(i)
{
  segments(x0 = path.2$x[i], y0 = path.2$y[i], 
           x1 = path.xy.2$x[i], y1 = path.xy.2$y[i],
           lwd = 1, lty = "dotted", col = "green")
} )

# 6. Eintragen der Spur in der xy-Ebene, also Verbindung aufeinanderfolgender Punkte in path.xy
lapply( X = (2:lgth), FUN = function(i)
{
  segments(x0 = path.xy.2$x[i-1], y0 = path.xy.2$y[i-1], 
           x1 = path.xy.2$x[i], y1 = path.xy.2$y[i],
           lwd = 1, lty = 1, col = "blue")
} )

# 7. Raumdiagonale berechnen:
s.max <- min(max(path[ , 1]), max(path[ , 2]), max(path[ , 3]))
s <- (1:s.max)
diagnl <- as.data.frame( cbind( s, s, s ) )

# 8. Raumdiagonale eintragen (schwarze Punkte)
scp$points3d(x = diagnl,
             type = "p",
             lwd = 1, col = "black")
             
# 9. Spur der Raumdiagonale eintragen (schwarz):
# Endpunkt der Raumdiagonale berechnen und mit Ursprung verbinden 
endpoint.2 <- scp$xyz.convert(x = diagnl[s.max, 1], y = diagnl[s.max, 2], z = 0)
segments(x0 = path.xy.2$x[1], y0 = path.xy.2$y[1], 
         x1 = endpoint.2$x, y1 = endpoint.2$y,
         lwd = 1, lty = 1, col = "black")