27 August 2009

Contourplot mit R

Neulich hab ich ja schon angedeutet, daß ich mal versuchen wollte, Höhenlinien mit dem Contourplot von R zu bauen. Gesagt, getan.
(Es folgt eigentlich weniger ein genaues "Rezept" als eine kurze Notiz - an mich? - ohne Gewähr):
Aus den srtm-Daten hab ich mir mit extractdata.pl einen Bereich um Trier ausgeschnitten, aus diesem Datensatz (trier.data) die Kopfzeilen per Hand gelöscht und diesen dann mit mosel <- as.matrix(read.table("trier.data", header=FALSE)) in R eingelesen. Um die Matrix zu "norden" (damit der plot dann auch "richtig rum" erscheint), muß diese gedreht werden (Funktionen dazu siehe hier): mosel = rotate270.matrix(mosel).
Danach geht's los mit dem Plot. Da hinter dem eigentlichen Contourplot die Höheninformationen als farbige Fläche angezeigt werden sollen*, bauen wir uns mit colpal = colorpanel(40,"green","red") zuerst eine Palette mit 40 Stufen aus den Farben grün bis rot. Mit image(mosel, col=colpal, axes=FALSE) kann man diesen "Hintergrund" dann schon mal plotten. (Die schöne Funktion filled.contour gibt es auch, diese erzeugt jedoch zwingend(?) eine Legende).
Den eigentlichen Contourplot legt man dann mit contour(mosel, nlevels=30, drawlabels = FALSE, axes = FALSE, add=TRUE) darüber. Ich hab hierbei mal auf alle Label und Achsenbeschriftungen verzichtet.
Den Plot kann man mit R auch gleich in ein Bitmap (z.B. mit jpeg) oder besser in ein SVG schreiben lassen: Zuerst ein device öffnen svg("mosel.svg"), die obigen Plots ausführen, danach das device schließen dev.off() und fertig ist das SVG.
Auf dem Bild sieht man das Ergebnis (SVG (~15MB) - mit GIMP nach Bitmap umgewandelt) für einen Ausschnitt (49,6522/6,5416 bis 49,8354/6,785), der als Matix aus srtm-Daten 292 Zeilen und 220 Spalten (also eigentlich wenige) umfasst. Als kleinen Test hab ich auch mal einen Contourplot einer kompletten srtm-Kachel (~140MB) erstellen lassen - wie man sich denken kann ist der Rechner dann bei 'vielen' Höhenstufen ganz gut ausgelastet, der import solcher Matrizen in R funktioniert aber noch ohne weiteres. Nur möchte ich nicht wissen, wie das Ganze aussieht, wenn man Contourplots von großen Gebieten erstellen möchte für die man Höheninformationen mit kleiner Rasterweite hat.

* um es mit Ivar Combrinck zu sagen: bunt "gibt die Anweisungen".

Labels: , , , , ,