REM Inizializzazione

g.mapset m=PERMANENT
r.mask rast=uds_cod1

REM importo e rasterizzo il layer litologia
v.in.ogr -e input="D:\GeoDati\lito\lito.shp" output="lito" --o
g.region rast=uds_cod1 res=5
v.clip -r input=lito output=litoc --o
v.to.rast input=litoc output=lito use=attr attribute_column=codnum --o

REM Creo un dato raster con il Corine Land Cover con tipi forestali al 4/5 livello
v.to.rast input=clc@PERMANENT output=clc_cod use=attr attribute_column=uso_cod --o

REM Creo un nuovo mapset diversita

g.mapset -c mapset=Diversita
g.region rast=uds_cod1 res=100

REM Biodiversità forestale
REM seleziono solo le celle con copertura forestale
r.mapcalc e="forest=if(clc_cod1==3,clc_cod,null())" --o
r.mask uds_cod1 maskcats=3

REM calcolo gli indici diversity ed interspersion in una finestra di 1000 metri pari a 11 celle
r.neighbors i=forest o=divfor1km method=diversity size=11 --o -c
r.neighbors i=forest o=intfor1km method=interspersion size=11 --o -c

REM calcolo indici morfologici

r.param.scale i=DEM o=geomorfo m=feature size=11 c=0.0002 --o

REM Calcolo Ecotopi

REM Analisi litologia
r.report m=lito units=k
REM Riclasificazione litologia
r.mapcalc e="c_lito=if(lito<8,1,if(lito<11,2,if(lito<15,3,4)))" --o

REM Analisi morfologia
r.report m=morfo units=k
REM Riclassificazione litologia
r.mapcalc "c_morfo=if(morfo<=2,1,if(morfo==3,2,3))"

REM Fisiotopo
r.cross i=c_lito,c_morfo o=fisiotopo -z --o

REM Ecotopi con Tipologie forestali
r.cross i=forest,c_morfo,c_lito o=ecotopi -z --o

REM Generalizzazione
r.neighbors i=ecotopi o=ecotopig size=5 method=mode -c --o

REM analisi diversità.
r.neighbors i=ecotopig o=div_ecotop size=11 method=diversity -c --o
r.neighbors i=ecotopig o=int_ecotop size=11 method=interspersion -c --o

REM La mappa è più comprensibile se classifichiamo la biodiversità con quantificatori linguistici 
REM p.es. ALTO MEDIO BASSO
REM Un possibile metodo per definire quantificatori relativi è basarsi sulla distribuzione di frequenza
REM il comando successivo calcola i 3 quantili
r.quantile i=div_ecotop q=3
r.mapcalc "div_ecotop_class=if(div_ecotop<4,4,if(div_ecotop<6,2,3))"
r.quantile i=int_ecotop q=5
r.mapcalc "int_ecotop_class=if(int_ecotop<32,1,if(int_ecotop<52,2,if(int_ecotop<65,3,if(int_ecotop<79,4,5))))"