[Fossgis-talk] [Grass-de] Raster-rechnen
Martin Schweizer
schweizermartin at students.unibe.ch
Mo Okt 12 18:18:13 CEST 2009
Hallo Markus und Liste
Vielen Dank für den Hinweis auf die Tutorials. Da war tatsächlich eine
Formel drin die ich gut gebrauchen kann.
(http://grass.osgeo.org/gdp/raster/mapcalc-algebra.pdf auf Seite 12)
Nun habe ich aber diese Formel eingegeben, es wird gerechnet, und am
Schluss habe ich eine Rasterkarte, bei welcher alle Zellen NULL-Werte
haben. Ich habe nun schon lange gesucht, aber einfach keinen Fehler
gefunden, deshalb hier mal die Frage: Findet jemand anders den Fehler?
Die Formel:
landw = sied+eval(x=eig+sied,\
if(x>(y=eig[-1,0]+sied[-1,0]),\
-.15*if(eig>y,sied,x-y),\
.15*if(eig[-1,0]>x,sied[-1,0],y-x))+\
if(x>(y=eig[1,0]+sied[1,0]),\
-.15*if(eig>y,sied,x-y),\
.15*if(eig[1,0]>x,sied[1,0],y-x))+\
if(x>(y=eig[0,-1]+sied[0,-1]),\
-.15*if(eig>y,sied,x-y),\
.15*if(eig[0,-1]>x,sied[0,-1],y-x))+\
if(x>(y=eig[0,1]+sied[0,1]),\
-.15*if(eig>y,sied,x-y),\
.15*if(eig[0,1]>x,sied[0,1],y-x))+\
if(x>(y=eig[-1,1]+sied[-1,1]),\
-.10*if(eig>y,sied,x-y),\
.10*if(eig[-1,1]>x,sied[-1,1],y-x))+\
if(x>(y=eig[1,1]+sied[1,1]),\
-.10*if(eig>y,sied,x-y),\
.10*if(eig[1,1]>x,sied[1,1],y-x))+\
if(x>(y=eig[1,-1]+sied[1,-1]),\
-.10*if(eig>y,sied,x-y),\
.10*if(eig[1,-1]>x,sied[1,-1],y-x))+\
if(x>(y=eig[-1,-1]+sied[-1,-1]),\
-.10*if(eig>y,sied,x-y),\
.10*if(eig[-1,-1]>x,sied[-1,-1],y-x)))
Und die Angaben zu den Karten:
GRASS 6.4.0RC5 (Landwirtschaft):~ > r.info map=eig
+----------------------------------------------------------------------------+
| Layer: eig Date: Mon Oct 12 18:07:15
2009
| Mapset: Ecthelias Login of Creator:
ecthelias
| Location:
Landwirtschaft
| DataBase:
/home/ecthelias/grassdata
| Title: ( eig
)
| Timestamp:
none
|----------------------------------------------------------------------------
|
| Type of Map: raster Number of Categories:
255
| Data Type:
FCELL
| Rows:
10417
| Columns:
10146
| Total Cells:
105690882
| Projection: UTM (zone
32)
| N: 5469466.69977 S: 4427807.08066979 Res:
99.99612356
| E: 1173129.00002557 W: 158549.86485219 Res:
99.99794354
| Range of data: min = -4.000000 max =
100.000000
|
| Data
Description:
| generated by
r.proj
|
|
Comments:
| r.proj input="Eignung52" location="Romanum5" mapset="Ecthelias"
outp\
| ut="eig"
method="nearest"
|
+----------------------------------------------------------------------------+
GRASS 6.4.0RC5 (Landwirtschaft):~ > r.info map=sied
+----------------------------------------------------------------------------+
| Layer: sied Date: Mon Oct 12 16:40:34
2009
| Mapset: Ecthelias Login of Creator:
ecthelias
| Location:
Landwirtschaft
| DataBase:
/home/ecthelias/grassdata
| Title: Labels ( Siedlungsflaeche2
)
| Timestamp:
none
|----------------------------------------------------------------------------
|
| Type of Map: raster Number of Categories:
0
| Data Type:
CELL
| Rows:
10417
| Columns:
10146
| Total Cells:
105690882
| Projection: UTM (zone
32)
| N: 5469466.69977 S: 4427807.08066979 Res:
99.99612356
| E: 1173129.00002557 W: 158549.86485219 Res:
99.99794354
| Range of data: min = 0 max =
20000000
|
| Data
Source:
| Vector Map: Siedlungsflaeche in mapset
Ecthelias
| Original scale from vector map:
1:1
|
| Data
Description:
| generated by
v.to.rast
|
|
Comments:
| v.to.rast input="Siedlungsflaeche" output="Siedlungsflaeche2"
use="a\
| ttr" type="point,line,area" layer=1 column="Flaeche" value=1
rows=4096
|
+----------------------------------------------------------------------------+
GRASS 6.4.0RC5 (Landwirtschaft):~ > r.info map=landw
+----------------------------------------------------------------------------+
| Layer: landw Date: Mon Oct 12 17:11:14
2009
| Mapset: Ecthelias Login of Creator:
ecthelias
| Location:
Landwirtschaft
| DataBase:
/home/ecthelias/grassdata
| Title: ( landw
)
| Timestamp:
none
|----------------------------------------------------------------------------
|
| Type of Map: raster Number of Categories:
255
| Data Type:
DCELL
| Rows:
10417
| Columns:
10146
| Total Cells:
105690882
| Projection: UTM (zone
32)
| N: 5469466.69977 S: 4427807.08066979 Res:
99.99612356
| E: 1173129.00002557 W: 158549.86485219 Res:
99.99794354
| Range of data: min = nan max =
nan
|
| Data
Description:
| generated by
r.mapcalc
|
|
Comments:
| sied + eval(x = eig + sied, if(x > (y = eig[-1,0] +
sied[-1,0]),
| -0.15 * if(eig > y, sied, x - y), 0.15 * if(eig[-1,0] >
x,
| sied[-1,0], y - x)) + if(x > (y = eig[1,0] + sied[1,0]), -0.15
*
| if(eig > y, sied, x - y), 0.15 * if(eig[1,0] > x, sied[1,0], y -
x))
| + if(x > (y = eig[0,-1] + sied[0,-1]), -0.15 * if(eig > y, sied, x
-
| y), 0.15 * if(eig[0,-1] > x, sied[0,-1], y - x)) + if(x > (y
=
| eig[0,1] + sied[0,1]), -0.15 * if(eig > y, sied, x - y), 0.15
*
| if(eig[0,1] > x, sied[0,1], y - x)) + if(x > (y = eig[-1,1]
+
| sied[-1,1]), -0.1 * if(eig > y, sied, x - y), 0.1 * if(eig[-1,1]
>
| x, sied[-1,1], y - x)) + if(x > (y = eig[1,1] + sied[1,1]), -0.1
*
| if(eig > y, sied, x - y), 0.1 * if(eig[1,1] > x, sied[1,1], y -
x))
| + if(x > (y = eig[1,-1] + sied[1,-1]), -0.1 * if(eig > y, sied, x
-
| y), 0.1 * if(eig[1,-1] > x, sied[1,-1], y - x)) + if(x > (y
=
| eig[-1,-1] + sied[-1,-1]), -0.1 * if(eig > y, sied, x - y), 0.1
*
| if(eig[-1,-1] > x, sied[-1,-1], y -
x)))
|
|
+----------------------------------------------------------------------------+
GRASS 6.4.0RC5 (Landwirtschaft):~ > g.region -p
projection: 1 (UTM)
zone: 32
datum: wgs84
ellipsoid: wgs84
north: 5469466.69977
south: 4427807.08066979
west: 158549.86485219
east: 1173129.00002557
nsres: 99.99612356
ewres: 99.99794354
rows: 10417
cols: 10146
cells: 105690882
GRASS 6.4.0RC5 (Landwirtschaft):~ >
Vielen Dank für jede Hilfe!
Viele Grüsse
Martin
Markus Neteler schrieb:
> 2009/10/5 Martin Schweizer <schweizermartin at students.unibe.ch>:
>
>> Hallo Markus
>>
>> Markus Neteler schrieb:
>>
>>> Hallo Martin,
>>>
>>>
> ...
>
>>> Du sagt vorne
>>> test=
>>> aber benutzt test auch in der Formel - wie soll das gehen?
>>>
>>>
>> Nun, die Karte "test", welche erstellt wird wird halt auch wieder in der
>> Formel verwendet. Soll heissen: Wenn auf der Karte "test" die soeben
>> erstellt wird, die betreffende Zelle bereits den Wert 5 besitzt, dann
>> bitte weiter zur nächsten Zelle (Dann ist die if-Bedingung nicht
>> erfüllt....).
>> Geht das aus irgend einem Grund nicht?
>>
>
> Ja genau, geht nicht (meines Wissens).
>
> Das Konzept ist
> neu=Funktion(alt)
>
> In der r.mapcalc Seite ist ein Tutorial verlinkt
> (etwas aelter, aber im Prinzip gueltig), was ein
> paar Modelle beinhaltet, die in Deine Richtung
> gehen.
>
> Gruesse
> Markus
>
>