[GRASS-SVN] r62486 - grass-addons/grass7/raster/r.rock.stability
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Oct 30 03:00:28 PDT 2014
Author: kikapu
Date: 2014-10-30 03:00:28 -0700 (Thu, 30 Oct 2014)
New Revision: 62486
Modified:
grass-addons/grass7/raster/r.rock.stability/r.rock.stability.py
Log:
r.rock.stability: upgrade
Modified: grass-addons/grass7/raster/r.rock.stability/r.rock.stability.py
===================================================================
--- grass-addons/grass7/raster/r.rock.stability/r.rock.stability.py 2014-10-30 08:12:54 UTC (rev 62485)
+++ grass-addons/grass7/raster/r.rock.stability/r.rock.stability.py 2014-10-30 10:00:28 UTC (rev 62486)
@@ -57,6 +57,7 @@
#%option
#% key: prefix
#% type: string
+#% gisprompt: new,raster,raster
#% description: Prefix for output maps
#% required: yes
#%end
@@ -131,57 +132,72 @@
aspect = 'aspects_' ,
zfactor = '1.0' ,
min_slp_allowed = '0.0' ,
+ overwrite = True ,
quiet = True)
# Correggi aspect per avere: N=0 E=90 S=180 W=270
grass.mapcalc("aspects_1 = if (aspects_ > 90, 450 - aspects_, 90 - aspects_)" ,
+ overwrite = True ,
quiet = True)
# calcolo A (valore assoluto)che vale SCIVOL=abs(immersdelgiunto-aspect1); RIBAL=abs(immersdelgiunto-aspect1-180)
grass.mapcalc("Asci = abs($imme - aspects_1)",
imme = imme ,
+ overwrite = True ,
quiet = True)
grass.mapcalc("Arib = abs(Asci - 180)" ,
+ overwrite = True ,
quiet = True)
# calcolo F1
# la formula originale e' stata in parte modificata per tener conto dei valori di A compresi tra 270 e 360 sfavorevoli x la stabilita'
grass.mapcalc("F1sci1 = if(Asci > 270, 0.64 - 0.006*atan(0.1*(abs(Asci - 360) - 17)), 0)" ,
+ overwrite = True ,
quiet = True)
grass.mapcalc("F1sci2 = if(Asci < 270, 0.64 - 0.006*atan(0.1*(Asci - 17)), 0)" ,
+ overwrite = True ,
quiet = True)
grass.mapcalc("F1sci = F1sci1 + F1sci2" ,
+ overwrite = True ,
quiet = True)
grass.mapcalc("F1rib1 = if(Arib > 270, 0.64 - 0.006*atan(0.1*(abs(Arib - 360) - 17)), 0)" ,
+ overwrite = True ,
quiet = True)
grass.mapcalc("F1rib2 = if(Arib < 270, 0.64 - 0.006*atan(0.1*(Arib - 17)), 0)" ,
+ overwrite = True ,
quiet = True)
grass.mapcalc("F1rib = F1rib1 + F1rib2" ,
+ overwrite = True ,
quiet = True)
# calcolo F2 (per il toppling vale sempre 1)
grass.mapcalc("F2sci = 0.56 + 0.005*(atan(0.17*abs($incl) - 5))" ,
incl = incl ,
+ overwrite = True ,
quiet = True)
# calcolo F3
grass.mapcalc("F3sci = -30 + 0.33*atan($incl - slopes_)" ,
incl = incl ,
+ overwrite = True ,
quiet = True)
grass.mapcalc("F3rib = -13 - 0.143*atan($incl + slopes_ - 120)" ,
incl = incl ,
+ overwrite = True ,
quiet = True)
# moltiplica F1xF2XF3
grass.mapcalc("Isciv = F1sci * F2sci * F3sci" ,
+ overwrite = True ,
quiet = True)
grass.mapcalc("Irib = F1rib * F3rib * 1" ,
+ overwrite = True ,
quiet = True)
# Calcola l'indice SMR
@@ -190,12 +206,14 @@
SMRsciv = prefix + '_SMRsciv' ,
RMR = RMR ,
F4 = f4 ,
+ overwrite = True ,
quiet = True)
grass.mapcalc("$SMRrib = Irib + $RMR + $F4" ,
SMRrib = prefix + '_SMRrib' ,
RMR = RMR ,
F4 = f4 ,
+ overwrite = True ,
quiet = True)
########################################################
############################SMR_wedge###################
@@ -213,9 +231,11 @@
imme2 = imme2 ,
incl = incl ,
incl2 = incl2 ,
+ overwrite = True ,
quiet = True)
grass.mapcalc("T1 = atan(T)" ,
+ overwrite = True ,
quiet = True)
#calcola plunge del cuneo P2 (sarebbe inclinazione)
grass.mapcalc("P=((sin($incl)*sin($imme)*sin($imme2-$imme))/(sin($incl)*cos($imme)*cos($incl2)-sin($incl2)*cos($imme2)*cos($incl)))*sin(T1)" ,
@@ -223,15 +243,20 @@
imme2 = imme2 ,
incl = incl ,
incl2 = incl2 ,
+ overwrite = True ,
quiet = True)
grass.mapcalc("P1=atan(P)" ,
+ overwrite = True ,
quiet = True)
grass.mapcalc("P2=abs(P1)" ,
+ overwrite = True ,
quiet = True)
#calcola trend del cuneo corretto (TB) sarebbe immersione
grass.mapcalc("TA=if(P1<0,T1+180,T1)" ,
+ overwrite = True ,
quiet = True)
grass.mapcalc("TB=if(TA<0,TA+360,TA)" ,
+ overwrite = True ,
quiet = True)
#calcola slope e aspect
#...calcolo A (valore assoluto)che vale SCIVOL=abs(immersdelgiunto-aspect1); RIBAL=abs(immersdelgiunto-aspect1-180)
@@ -244,9 +269,11 @@
output='F1sci_',
rules='-',
stdin=reclass_rules,
+ overwrite = True ,
quiet = True)
# creo una mappa con valore dell'inclinazione
grass.mapcalc('mappa=P2+slopes_*0' ,
+ overwrite = True ,
quiet = True)
# calcolo di F2
reclass_rules2 = "0 thru 20 = 0.15\n20.1 thru 30 = 0.4\n30.1 thru 35 = 0.7\n35.1 thru 45 = 0.85\n45.1 thru 90 = 1"
@@ -255,10 +282,12 @@
output='F2',
rules='-',
stdin=reclass_rules2,
+ overwrite = True ,
quiet = True )
# Calcolo di F3
# Calcolo C che vale SCIVOL=inclgiunto - slopes;
grass.mapcalc("Csci=P2-slopes_" ,
+ overwrite = True ,
quiet = True)
# Reclass per F3 scivolamento
reclass_rules3 = "180 thru 9 = 0\n10 thru 0.1 = -6\n0 = -25\n-0.1 thru -9 = -50\n-10 thru -180 = -60"
@@ -266,6 +295,7 @@
input='Csci' ,
output='F3sci_' ,
rules='-' ,
+ overwrite = True ,
stdin=reclass_rules3)
grass.mapcalc('Isciv_=F1sci_*F3sci_*F2' ,
quiet = True)
@@ -273,11 +303,10 @@
SMRsciv = prefix + '_SMRscivWedge' ,
RMR = RMR ,
F4 = f4 ,
+ overwrite = True ,
quiet = True)
grass.run_command('g.remove' ,
- flags='f',
- type='rast',
- name=(
+ rast=(
'Asci' ,
'Arib' ,
'F1sci_' ,
@@ -306,70 +335,87 @@
# calcolo DELTA=aspect-immersione
grass.mapcalc("delta = aspects_1 - $imme",
imme = imme ,
+ overwrite = True ,
quiet = True)
# calcola A=cosenoDELTA
grass.mapcalc("A = cos(delta)" ,
+ overwrite = True ,
quiet = True)
# calcola B=tan(inclinazione)
grass.mapcalc("B = tan($incl)" ,
incl = incl ,
+ overwrite = True ,
quiet = True)
# calcola C=A*B
grass.mapcalc("C = A * B" ,
+ overwrite = True ,
quiet = True)
# calcola AP
grass.mapcalc("AP = atan(C)" ,
+ overwrite = True ,
quiet = True)
# additional condition per scivolamento (dip del pendio deve essere < di AP+5), AP deve essere < di 85
grass.mapcalc("AP1 = slopes_ - AP - 5" ,
+ overwrite = True ,
quiet = True)
grass.mapcalc("consci1 = if(AP1 > 0, 0, null())" ,
+ overwrite = True ,
quiet = True)
grass.mapcalc("consci2 = if(AP < 85, 0, null())" ,
+ overwrite = True ,
quiet = True)
# additional condition per ribaltamento (AP deve essere> di 85)
grass.mapcalc("conrib = if(AP > (-85), 0, null())" ,
+ overwrite = True ,
quiet = True)
# ANALISI SCIVOLAMENTO PLANARE AP * 0.0113
grass.mapcalc("membro_sci = 0.0113 * AP" ,
+ overwrite = True ,
quiet = True)
grass.mapcalc("sci = membro_sci - $TC" ,
TC = TC ,
+ overwrite = True ,
quiet = True)
# se il risultato ottenuto e' minore di zero non ho cinematismo, maggiore e' il valore piu' alta sara' la propensione al distacco
grass.mapcalc("cinem_sci = if(sci >= 0, 1, null())" ,
TC = TC ,
+ overwrite = True ,
quiet = True)
grass.mapcalc("$SSPCsciv = cinem_sci * sci + consci1 + consci2" ,
SSPCsciv = prefix + '_SSPCsciv' ,
+ overwrite = True ,
quiet = True)
# ANALISI RIBALTAMENTO 0.0087*(-90-AP+INCLINAZIONE)
grass.mapcalc("membro_RIB = 0.0087 * ((-90) - AP + $incl)" ,
incl = incl ,
+ overwrite = True ,
quiet = True)
grass.mapcalc("rib = membro_RIB - $TC" ,
TC = TC ,
+ overwrite = True ,
quiet = True)
# se il risultato ottenuto e' minore di zero non ho cinematismo, maggiore e' il valore piu' alta sara' la propensione al distacco
grass.mapcalc("cinem_rib = if(rib >= 0, 1, null())" ,
+ overwrite = True ,
quiet = True)
grass.mapcalc("$SSPCrib = cinem_rib * rib + conrib" ,
SSPCrib = prefix + '_SSPCrib' ,
+ overwrite = True ,
quiet = True)
# processo per sostituire i valori di cella nulli con 0 (passaggio importante per trovare poi il minimo
grass.run_command("r.null" ,
map = prefix + '_SSPCsciv' ,
null = 0 ,
+ overwrite = True ,
quiet = True)
grass.run_command("r.null" ,
map = prefix + '_SSPCrib' ,
null = 0 ,
+ overwrite = True ,
quiet = True)
#elimino mappe
grass.run_command('g.remove' ,
- flags='f',
- type='rast',
- name=(
+ rast=(
'slopes_' ,
'aspects_' ,
'aspects_1' ,
@@ -404,9 +450,7 @@
quiet = True)
else:
grass.run_command('g.remove' ,
- flags='f',
- type='rast',
- name=(
+ rast=(
'slopes_' ,
'aspects_' ,
'aspects_1' ,
More information about the grass-commit
mailing list