Kartenalgebra mit r.mapcalc

Figure 16: "`Moving-window"'-Methode bei Kartenalgebra im Rasterformat: Der Rasterzellenwert X in der neuen Karte berechnet sich hier aus einer 3x3-Matrix der alten Karte.

Kartenverschneidungen auf Rasterebene bieten eine enorme Bandbreite an Analysemöglichkeiten. In GRASS kann mit dem r.mapcalc-Modul Kartenalgebra durchgeführt werden. Es arbeitet auch mit der "`moving window"'-Technologie, dabei bewegt sich ein Fenster definierbarer Größe sukzessive durch alle Spalten und Reihen und verrechnet zellenweise Rasterkarten miteinander (vgl. Abb. 16). Ab GRASS 5.0.x wird zwischen NULL ("`no data"', also "`nichts"') und 0 (ZERO) unterschieden.

Prinzipiell wird in $ r.mapcalc auf folgende Art gerechnet:

   neuekarte = Wert/karte1 $<$formel$>$ [karte2, ...]

Links vor dem Gleichzeichen steht der Name der zu erzeugenden Rasterkarte, rechts die Berechnungsformel. Das Modul kann also "`intuitiv"' benutzt werden.

Es gibt folgende Operatoren in r.mapcalc:

 
Folgende Funktionen stehen in r.mapcalc zur Verfügung:

 

Zusätzlich gibt es einige interne Variablen in r.mapcalc, die Werte beziehen sich auf das "`moving window"', das die Berechnung durchführt:

 

Der Wert "`keine Daten"' (NULL, "`no data"') wird mit "`null()"' angegeben. NULL ist von 0 (ZERO) verschieden.

Einfache Beispiele für die Grundrechenarten: Auf einfache Art und Weise können die geographisch korrespondierenden Zellenwerte von Karten miteinander addiert werden.

Beispiel 1: Addiere Gebäudestrukturen (als Gebäudehöhen in Karte "`gebaeude"' gespeichert) zu einem Höhenmodell, so dass die Gebäude im Gelände "`stehen"':

   dgm mit haeusern = gebaeude + hoehenmodell

Beispiel 2: Berechung eines gewichteten Mittels aus zwei Karten (Dezimalpunkte, damit "`neue"-karte"' als Fließkommakarte erzeugt wird):

   neuekarte = (5. $\ast$ karte1 + 3. $\ast$ karte2)/8.2

Einfache Beispiele für if-Bedingungen: ("`karte"' ist ein Kartenname, alle anderen Zeichen können Kartennamen oder Werte sein, die Rechnung erfolgt pixelweise):

1. if karte = a then b else c

   in r.mapcalc: neuekarte = if((karte == a),b,c)

2. if karte $<$$>$ a then b else c

   in r.mapcalc: neuekarte = if((karte != a),b,c)

3. if karte $>$= a then b else c

   in r.mapcalc: neuekarte = if((karte $>$= a),b,c)

4. if a $<$= karte $<$= b then c else d

   in r.mapcalc: neuekarte = if(((karte $>$= a) && (karte $<$= b)),c,d)

Hier nun weitere Beispiele für komplexere Datenselektionen in Rasterkarten:

a) Selektiere aus "`karte"' die Werte 1 und 2 und speichere in neuer Karte unter Beibehaltung der selektierten Werte, 0 sonst:

   neuekarte = if((karte==1),1,0) + if((karte==2),2,0)

b) Selektiere aus "`karte"' die Werte 1 und 2 und speichere in neuer Karte als Binärkarte ab (nur Werte 0 und 1 (also "`ja"' und "`nein"')):

   neuekarte = if((karte==1),1,0) || if((karte==2),2,0)

c) Berücksichtigung von NULL-Werten ("`no data"', nichts):

 i. Umwandlung ausschließlich von NULL-Werten in einen neuen Wert:

if karte=NULL then 3 else karte (wenn "`no data"', dann 3 sonst Pixelwert übernehmen)

   in r.mapcalc: neuekarte = if(isnull(karte),3,karte)

 ii. Umwandlung aller Zellenwerte kleiner 5 in der Karte in den NULL-Wert, alle größeren in 5:

   in r.mapcalc: neuekarte = if(karte$<$5,null(),5)

Eine weitere interessante Möglichkeit bietet die Angabe relativer Koordinaten, die für das "`wandernde Fenster"' gültig sind. Damit können auch Nachbarzellen mit in die Rechnung integriert bzw. ein größeres/unsymmetrisches "`moving window"' als die übliche 3x3-Matrix definiert werden. Beispiel:

   neuekarte = altekarte[1,1] + altekarte[-1,-1]

benutzt nur die Zellen oben rechts [1,1] und unten links [-1,-1] in einer 3x3-Matrix zur Berechnung der neuen Karte. Diese Möglichkeit lässt sich auch bei verschiedenen Eingangskarten einsetzen. Mit "`row()"' und "`col()"' können die aktuellen Zeilen- und Spaltenwerte in der Formel integriert werden, mit "`x()"' und "`y()"' die laufenden Koordinaten des "`moving window"'.

Zur einfacheren Editierung (Verwendung der Cursortasten) kann r.mapcalc auch kommandozeilenorientiert benutzt werden. Dabei ist die Formel in Anführungsstriche zu setzen. Ein Beispiel für die Erzeugung eines Bildmosaiks: wenn Koordinaten in x- und y-Richtung kleiner als angegebener Wert, dann LANDSAT-Kanal 1 (tm1) ausgeben, sonst an diesem Pixel die Topographische Karte (tk25) darstellen:

$ r.mapcalc 'neuekarte=if((x()$<$3571000 && y()$<$5767000),tm1,tk25)'

Hier wurde von Variablen Gebrauch gemacht, die die laufenden Koordinaten des "`moving windows"' speichern. Problematisch ist hier nur, dass die Farben der Originalkarten nicht erhalten bleiben. Auf ähnliche Weise kann auch schnell ein Rechteck mit definierten Randkoordinaten und bestimmten Zellenwerten erzeugt werden.

Eine praktische Funktion ist "`eval()"'. Damit lassen sich Zwischenwerte speichern, ohne neue Karten erzeugen zu müssen und sich so mehrere Schritte in eine Formel bringen. Die einzelnen Formeln werden innerhalb der "`eval()"'-Umgebung mit Komma getrennt. Die letzte Anweisung wird hier immer als Ergebnis gespeichert. Interessant ist diese Funktion beispielsweise bei Abfragen von Fließkommakarten, wenn nicht exakte Bereiche definiert werden sollen. In diesem Beispiel werden die Kartenwerte erst gerundet, in "`a"' zwischengespeichert und dann über eine Bedingung abgefragt:

(wenn gerundete Zellenwerte der Karte =3, dann 3 erhalten, sonst NULL setzen)

neuekarte = eval( a=round(karte), if(a==3,3,null()) )

Als komplexeres Beispiel soll noch einmal eine Bereichsselektion in einer Fließkommakarte gezeigt werden, bei der die NULL-Werte ("`no data"') erhalten bleiben, Speicherung der Zwischenwerte in die Variablen "`zw1"' und "`zw2"':

if a $<$= karte $<$= b then c else d mit Erhalt von NULL (in einer Befehlszeile eingeben):

neuekarte = eval (zw1=round(karte), zw2=if(((zw1 $>$= a) && (zw2 $<$= b)),c,d),

  if(isnull(karte), karte, zw2))

Zur Veranschaulichung noch einmal als Zahlenbeispiel:

if 4 $<$= round(karte) $<$= 5 then 20 else 66 mit Erhalt von NULL (in einer Befehlszeile eingeben):

neuekarte = eval (zw1=round(karte), zw2=if(((zw1 $>$= 4) && (zw1 $<$= 5)),20,66),

  if(isnull(karte),karte,zw2))

Es muss selbstverständlich nicht vor Berechnungen bei Fließkommakarten gerundet werden, Sie können auch Bereiche oder explizite Werte angeben.

Beispiel: Erzeugung einer in Nordwest-Richtung geneigten Ebene mit Starthöhe 100m:

neuekarte = 100 + row() + col()

Weitere Beispiele zu r.mapcalc finden sich bei ALBRECHT 1992 und dem Tutorial zu r.mapcalc von SHAPIRO und WESTERVELT 1992.