[GRASS-SVN] r31204 - in grass-addons/raster: . r.inund.fluv
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri May 2 08:24:52 EDT 2008
Author: robertomarzocchi
Date: 2008-05-02 08:24:52 -0400 (Fri, 02 May 2008)
New Revision: 31204
Added:
grass-addons/raster/r.inund.fluv/
grass-addons/raster/r.inund.fluv/2d_path.f90
grass-addons/raster/r.inund.fluv/README.txt
grass-addons/raster/r.inund.fluv/clean_inundation.f90
grass-addons/raster/r.inund.fluv/correction_from_path.f90
grass-addons/raster/r.inund.fluv/find_main_channel.f90
grass-addons/raster/r.inund.fluv/install.sh
grass-addons/raster/r.inund.fluv/r.inund.fluv
grass-addons/raster/r.inund.fluv/r.inund.fluv.html
grass-addons/raster/r.inund.fluv/r.inund.fluv_62
Log:
Added: grass-addons/raster/r.inund.fluv/2d_path.f90
===================================================================
--- grass-addons/raster/r.inund.fluv/2d_path.f90 (rev 0)
+++ grass-addons/raster/r.inund.fluv/2d_path.f90 2008-05-02 12:24:52 UTC (rev 31204)
@@ -0,0 +1,1182 @@
+module dd
+integer,parameter::ferma_bisezione=15
+integer::n_no_null
+real, allocatable, dimension(:)::val, iind, jind
+real::output
+end module dd
+
+module dati
+!****************************************************************************************************************************************
+!real, parameter::=2
+!il programma è pensato per una risoluzione del file raster 2x2 se la risoluzione fosse diversa MODIFICARE il parametro res
+!****************************************************************************************************************************************
+real::nord, sud, est, ovest
+integer:: righe, colonne, i, j, righe2, colonne2
+integer, allocatable, dimension(:,:):: bordi_alveo
+real, allocatable, dimension(:,:):: quota_terreno, reticolo2d, inondazioneT
+real:: Ebordo, Nbordo, EEE,NNN, zzz, res, res2
+end module dati
+
+
+
+
+!****************************************************************************************************************************************
+!il programma calcola le direzioni privilegiate dell'acqua al di fuori dell'alveo ottenuto con il precedente modulo trova_alveo.f90
+! a tale scopo utilizza in input: - il file contenente i limiti della regione
+! - il file contenente il profilo
+! - il file limiti_alveo_vect
+! - il DTM
+!fornisce in output per il profilo considerato una carta contenente le direzioni privilegiate dell'acqua in base alla massima pendenza
+!****************************************************************************************************************************************
+
+
+program reticolo_2d_fuori_alveo
+Use dati
+Use dd
+implicit none
+
+integer:: ii, f, iend, nrighe, h1, h2, h3, h4, pixel_temp_i, pixel_temp_j
+integer:: nsez, nbordi, nscelta, cont
+!integer, allocatable, dimension(:,:):: bordi_alveo
+!real, allocatable, dimension(:,:):: quota_terreno
+real(kind=8)::E1, E1dopo, N1, N1dopo, quotal, temp, temp1, temp2, z, temp_x, temp_y, E_temp, N_temp, dist
+real(kind=8):: num, den, x, y
+real(kind=8),dimension(2):: El_ds, El_sx, Nl_ds, Nl_sx
+real, allocatable, dimension(:)::profilo, E, N
+character(len=40)::xx, tab, tab1,vv
+character(40)::nomefile1,nomefile2,nomefile3,nomefile4
+character(40)::nomefile5,nomefile6,nomefile7,nomefile8,nomefile9
+open(unit=40, status='old', file='nomefile.txt')
+read(40,*)nomefile1,nomefile2,nomefile3,nomefile4,nomefile5,&
+nomefile6,nomefile7,nomefile8,nomefile9
+close(40)
+
+!50 region20
+!51 profilo
+!52 dtm20x20'
+!53 limiti_alveo_vect
+!54 inondazione_step3
+!55 region2.txt
+! dtm2x2 nomefile8 (vd dopo)
+!60 bordi_alveo
+!61 reticolo
+
+
+
+open(unit=50, status='old', file=nomefile4)
+open(unit=51, status='old', file=nomefile1)
+open(unit=52, status='old', file=nomefile5)
+open(unit=53, status='old', file=nomefile6)
+open(unit=54, status='old', file=nomefile2)
+open(unit=55, status='old', file=nomefile7)
+!open(unit=56, status='old', file='grass_script/dtm2x2')
+open(unit=60, file=nomefile9)
+open(unit=61, file=nomefile3)
+
+
+!leggo le dimensioni della regione utili per allocare le matrici in cui memorizzare carta di input e di output
+do i=1,13
+ read(50,*)xx
+ if (xx=='north:') then
+ backspace(50)
+ read(50,*)xx, nord
+ !write(60,'(a6,i)') trim(xx), nord
+ !write(61,'(a6,i)') trim(xx), nord
+ ! write(60,'(a6,f10.2)') trim(xx), nord
+ write(61,'(a6,f10.2)') trim(xx), nord
+ !write(60,*) trim(xx), nord
+ !write(61,*) trim(xx), nord
+ else if (xx=='nsres:') then
+ backspace(50)
+ read(50,*)xx, res
+ else if (xx=='south:') then
+ backspace(50)
+ read(50,*)xx, sud
+ !write(60,'(a6,f10.2)')trim(xx), sud
+ write(61,'(a6,f10.2)')trim(xx), sud
+ !write(60,*)trim(xx), sud
+ !write(61,*)trim(xx), sud
+ else if (xx=='east:') then
+ backspace(50)
+ read(50,*)xx, est
+ !write(60,'(a5,f10.2)') trim(xx), est
+ write(61,'(a5,f10.2)') trim(xx), est
+ !write(60,*)trim(xx), est
+ !write(61,*)trim(xx), est
+ else if (xx=='west:') then
+ backspace(50)
+ read(50,*)xx, ovest
+ !write(60,'(a5,f10.2)')trim(xx), ovest
+ write(61,'(a5,f10.2)')trim(xx), ovest
+ !write(60,*)trim(xx), ovest
+ !write(61,*)trim(xx), ovest
+ else if (xx=='rows:') then
+ backspace(50)
+ read(50,*)xx, righe
+ !write(60,'(a5,i5)')trim(xx), righe
+ write(61,'(a5,i5)')trim(xx), righe
+ !write(60,*)trim(xx), righe
+ !write(61,*)trim(xx), righe
+ else if (xx=='cols:') then
+ backspace(50)
+ read(50,*)xx, colonne
+ !write(60,'(a5,i5)')trim(xx), colonne
+ write(61,'(a5,i5)')trim(xx), colonne
+ !write(60,*)trim(xx), colonne
+ !write(61,*)trim(xx), colonne
+ endif
+enddo
+write(6,*) nord, sud, est, ovest, righe, colonne,res
+rewind(50)
+close(50)
+
+allocate(quota_terreno(righe,colonne))
+allocate(bordi_alveo(righe,colonne))
+allocate(reticolo2d(righe,colonne))
+
+!leggo il DTM
+read(52,*) ((quota_terreno(i,j), j=1,colonne), i=1,righe)
+write(6,*) "Read the Digital Terrain Model lower resolution"
+rewind(52)
+close(52)
+
+!scrivo 0 sulla matrice bordi_alveo e reticolo2d
+do j=1,colonne
+ do i=1,righe
+ bordi_alveo(i,j)=0
+ reticolo2d(i,j)=0
+ enddo
+enddo
+
+
+!guardo quant'è lungo il file che contiene il profilo
+do i=1,10000000
+ read(51,*,iostat=iend)xx
+ if(iend /= 0)exit
+enddo
+if(iend > 0)then
+ write(6,*)'Error reading the file with water surface profile :',iend
+ stop
+endif
+nrighe=i-1
+write(6,*) 'The file with water surface profile has',nrighe,' rows'
+rewind(51)
+
+
+!guardo quant'è lungo il file che contiene i bordi sezioni
+do i=1,10000000
+ read(53,*,iostat=iend)xx
+ if(iend /= 0)exit
+enddo
+if(iend > 0)then
+ write(6,*)"Error reading the file with river's boundary:",iend
+ stop
+endif
+nbordi=i-1
+write(6,*) "The file with river's boundary has ",nbordi,'rows'
+rewind(53)
+
+
+!leggo le dimensioni della regione utili per allocare le matrici in cui memorizzare carta di input e di output
+do i=1,13
+ read(55,*)vv
+ if (vv=='nsres:') then
+ backspace(55)
+ read(55,*)vv, res2
+ else if (vv=='rows:') then
+ backspace(55)
+ read(55,*)vv, righe2
+ else if (vv=='cols:') then
+ backspace(55)
+ read(55,*)vv, colonne2
+ endif
+enddo
+write(6,*) righe2, colonne2,res2
+rewind(55)
+close(55)
+
+
+allocate(inondazioneT(righe,colonne))
+!allocate(terreno(righe2,colonne2))
+
+!leggo il DTM 2x2
+!read(56,*) ((terreno(i,j), j=1,colonne2), i=1,righe2)
+!write(6,*) 'letto il DTM'
+!rewind(56)
+!close(56)
+
+call matrice_sparsa (nomefile8,righe2,colonne2)
+write(6,*) "Read the Digital Terrain Model with Compressed Sparse Row (CSR) algorithm"
+
+!leggo la carta inondazione
+read(54,*) ((inondazioneT(i,j), j=1,colonne), i=1,righe)
+write(6,*) 'Read the raster map of inondation of step 3'
+rewind(54)
+close(54)
+
+
+
+
+
+!leggo il file contenente il profilo e se sono nella mia regione traccio delle rette congiungenti i limiti alveo trovati in ogni sezione col precedente modulo
+f=0
+i=0
+cont=0
+1000 i=i+1
+if (i==nrighe) then
+ go to 1001
+endif
+read(51,*)E1, N1
+!write(6,*)i, E1, N1
+read(51,*)E1dopo, N1dopo
+!write(6,*)i, E1dopo, N1dopo
+backspace(51)
+if (E1<ovest .or. E1>est .or. N1>nord .or. N1<sud) then
+ go to 1000
+endif
+f=f+1
+backspace(51)
+read(51,*) E1, N1, z
+!write(6,*)cont, i, E1, N1, z
+if(cont==nbordi) then
+ go to 1000
+endif
+read(53,*) El_ds(1), Nl_ds(1), quotal, temp1, h1
+cont=cont+1
+if(cont==nbordi) then
+ go to 1000
+endif
+read(53,*) El_sx(1), Nl_sx(1), quotal, temp1, h2
+cont=cont+1
+if(cont==nbordi) then
+ go to 1000
+endif
+!1write(6,*)El_ds(1), Nl_ds(1), quotal, temp1, h1,h2
+!controllo che ci sia corrispondenza tra il limite alveo a destra e a sinistra (in realtà vale solo per il primo punto)
+if(h1/=h2)then
+ backspace(53)
+ cont=cont-1
+ go to 1000
+endif
+if(cont==nbordi) then
+ go to 1000
+endif
+123 read(53,*) El_ds(2), Nl_ds(2), quotal,temp,h3
+cont=cont+1
+if(cont==nbordi) then
+ go to 1000
+endif
+!write(6,*)El_ds(2), Nl_ds(2), quotal,temp,h3
+read(53,*) El_sx(2), Nl_sx(2), quotal,temp,h4
+cont=cont+1
+!write(6,*)cont, f, El_sx(2), Nl_sx(2), quotal,temp,h3,h4
+!controllo che ci sia corrispondenza tra il limite alveo a destra e a sinistra in caso contario "salto quel punto" con attenzione
+if(h3/=h4)then
+ backspace(53)
+ cont=cont-1
+ i=i+1
+ f=f+1
+ read(51,*) x
+ read(51,*) E1dopo, N1dopo
+ backspace(51)
+ if(cont==nbordi) then
+ go to 1000
+ endif
+ go to 123
+ !read(53,*) El_ds(2), Nl_ds(2), quotal,temp,h3
+ !cont=cont+1
+ !if(cont==nbordi) then
+ ! go to 1000
+ !endif
+ !read(53,*) El_sx(2), Nl_sx(2), quotal,temp,h4
+ !cont=cont+1
+ !write(6,*) f, El_sx(2), Nl_sx(2), quotal,temp,h3,h4
+endif
+backspace(53)
+cont=cont-1
+backspace(53)
+cont=cont-1
+
+
+!ora devo tracciare le varie rette, prima però devo controllare che effettivamente (El_ds(i),Nl_ds(i)) siano entrambe alla destra della centerline
+!per fare ciò controllo che il segmento che unisce (E1,N1) a (E1dopo,N1dopo) e quello che unisce i limiti alveo non si incrocino,
+!qualora si incrociassero inverto ds(2) e sx(2)
+num=N1-(N1dopo-N1)/(E1dopo-E1)*E1-Nl_ds(1)+(Nl_ds(2)-Nl_ds(1))/(El_ds(2)-El_ds(1))*El_ds(1)
+den=(Nl_ds(2)-Nl_ds(1))/(El_ds(2)-El_ds(1))-(N1dopo-N1)/(E1dopo-E1)
+x=num/den
+y=N1+(N1dopo-N1)/(E1dopo-E1)*(x-E1)
+!if(f==41)then
+! write(6,*)x,y,E1,N1,E1dopo,N1dopo,El_ds(1),Nl_ds(1),El_ds(2),Nl_ds(2)
+!endif
+
+!write(6,*) num, den, x, y
+!guardo se si incrociano
+!1/2
+if ( y<=Nl_ds(2) .and.y>=Nl_ds(1) .and. El_ds(1)<=x .and. x<=El_ds(2)) then
+ !write(6,*) 'interseca --> correggo1'
+ go to 2000
+!1\2
+else if ( El_ds(1)<=x .and. x<=El_ds(2) .and. Nl_ds(2)<=y .and. y<=Nl_ds(1)) then
+ !write(6,*) 'interseca --> correggo2'
+ go to 2000
+!2\1
+else if (El_ds(2)<=x .and. x<=El_ds(1) .and. y<=Nl_ds(2) .and. y>=Nl_ds(1)) then
+ !write(6,*)'interseca --> correggo3'
+ go to 2000
+!2/1
+else if (El_ds(2)<=x .and. x<=El_ds(1) .and. Nl_ds(2)<=y .and. y<=Nl_ds(1)) then
+ !write(6,*) 'interseca --> correggo4'
+ go to 2000
+else
+ !write(6,*)'bene'
+ go to 2001
+endif
+
+
+2000 temp1=Nl_ds(2)
+temp2=Nl_sx(2)
+Nl_ds(2)=temp2
+Nl_sx(2)=temp1
+temp1=El_ds(2)
+temp2=El_sx(2)
+El_ds(2)=temp2
+El_sx(2)=temp1
+go to 2001
+
+
+
+!a questo punto sono sicuro che i limiti ds e sx siano entrambi dallo stesso lato rispetto alla centerline e procedo col tracciare i bordi dell'alveo
+!interpolando linearmente fra 1 punto e il successivo
+!DESTRA
+2001 write (60,'(a6)')'L 2 1'
+write (60,*)El_ds(1),Nl_ds(1)
+write (60,*)El_ds(2),Nl_ds(2)
+write (60,*)' 1 1'
+dist=sqrt((El_ds(2)-El_ds(1))**2.+(Nl_ds(2)-Nl_ds(1))**2.)
+E_temp=El_ds(1)
+N_temp=Nl_ds(1)
+temp=0
+20011 x=res/50.
+temp=temp+x
+E_temp=E_temp+x*(El_ds(2)-El_ds(1))/dist
+N_temp=N_temp+x*(Nl_ds(2)-Nl_ds(1))/dist
+!write(6,*)E_temp,N_temp,dist
+if(temp>=dist)then
+ !write(6,*)'fin qua tutto ok'
+ go to 2002
+endif
+temp_x=int(E_temp)-ovest
+temp_y=nord-int(N_temp)
+pixel_temp_i=1+int(temp_y/res)
+pixel_temp_j=1+int(temp_x/res)
+!write(6,*) pixel_temp_i, pixel_temp_j
+bordi_alveo(pixel_temp_i,pixel_temp_j)=1
+go to 20011
+
+
+!SINISTRA
+2002 write (60,'(a6)')'L 2 1'
+write (60,*)El_sx(1),Nl_sx(1)
+write (60,*)El_sx(2),Nl_sx(2)
+write (60,*)' 1 2'
+dist=sqrt((El_sx(2)-El_sx(1))**2.+(Nl_sx(2)-Nl_sx(1))**2.)
+E_temp=El_sx(1)
+N_temp=Nl_sx(1)
+temp=0
+20021 x=res/50.
+temp=temp+x
+E_temp=E_temp+x*(El_sx(2)-El_sx(1))/dist
+N_temp=N_temp+x*(Nl_sx(2)-Nl_sx(1))/dist
+if(temp>=dist)then
+ go to 1000
+endif
+temp_x=int(E_temp)-ovest
+temp_y=nord-int(N_temp)
+pixel_temp_i=1+int(temp_y/res)
+pixel_temp_j=1+int(temp_x/res)
+bordi_alveo(pixel_temp_i,pixel_temp_j)=1
+go to 20021
+
+1001 rewind(51)
+rewind(53)
+
+
+!ora mi memorizzo il profilo in un'array opportunamente creato
+allocate(profilo(nrighe))
+allocate(E(nrighe))
+allocate(N(nrighe))
+
+
+do i=1,nrighe
+ profilo(i)=0
+enddo
+
+
+f=0
+i=0
+3000 i=i+1
+if (i==nrighe) then
+ go to 4000
+endif
+read(51,*)E1, N1, z
+!write(6,*)i, E1, N1
+!write(6,*)i, E1dopo, trim(tab), N1dopo
+if (E1<ovest .or. E1>est .or. N1>nord .or. N1<sud) then
+ go to 3000
+else
+ f=f+1
+ E(f)=E1
+ N(f)=N1
+ profilo(f)=z
+ go to 3000
+endif
+
+
+4000 nsez=f
+write(6,*) "There are", nsez, "section"
+rewind(51)
+
+
+do ii=1,nbordi
+ !write(6,*) ii, 'ok', nbordi
+ read(53,*) Ebordo, Nbordo, quotal, temp1, h1
+ ! write(6,*) Ebordo, Nbordo, quotal, temp1, h1
+ if(h1/=1 .or. h1/=nsez) then
+ EEE=E(h1)
+ NNN=N(h1)
+ zzz=profilo(h1)
+ !write(6,*)ii, 'fin qua ok'
+ call muoviti
+ endif
+enddo
+
+
+
+write(6,*) 'End of computational step → WRITING OUTPUT'
+!write(60,*) ((bordi_alveo(i,j), j=1,colonne), i=1,righe)
+write(61,*) ((reticolo2d(i,j), j=1,colonne), i=1,righe)
+
+write(6,*) 'END of Fortran Code → OK'
+end program reticolo_2d_fuori_alveo
+
+
+
+
+
+
+
+
+
+
+
+!****************************************************************************************************************************************
+!SUBROUTINE DIREZIONE
+!****************************************************************************************************************************************
+subroutine muoviti
+USE dati
+USE dd
+implicit none
+
+real::EST1,EST2,NORD1,NORD2,LIVELLO
+!real,intent(IN)::EST1
+!real,intent(IN)::EST2
+!real,intent(IN)::NORD1
+!real,intent(IN)::NORD2
+!real,intent(IN)::LIVELLO
+integer::iii, i_temp, j_temp, iiii, jjjj
+real::vettore_direzione, quota_min, temp, dist_temp, xxx, E_t, N_t, passo, xxxx
+character(1)::dir
+!(EEE,Ebordo,NNN,Nbordo,zzz)
+!(EST1,EST2,NORD1,NORD2,LIVELLO)
+
+EST1=EEE
+NORD1=NNN
+EST2=Ebordo
+NORD2=Nbordo
+LIVELLO=zzz
+
+!****************************************************************************************************************************************
+!SOLO IN ENTRATA
+quota_min=0
+i=1+int((nord-NORD2)/res)
+j=1+int((EST2-ovest)/res)
+vettore_direzione=((NORD2-NORD1)/(EST2-EST1))
+!assegno dir positiva qualora y(i+1)>y(i) viceversa dir negativa
+if(NORD2>=NORD1)then
+ dir='+'
+else
+ dir='-'
+endif
+dist_temp=sqrt((NORD2-NORD1)**2+(EST2-EST1)**2)
+iii=0
+1111 xxx=res
+temp=EST2+xxx*(EST2-EST1)/dist_temp
+j=1+int((temp-ovest)/res)
+temp=NORD2+xxx*(NORD2-NORD1)/dist_temp
+i=1+int((nord-temp)/res)
+!write(6,*)i, j, righe, colonne
+if(i<1 .or. j<1 .or. i>righe .or. j>colonne) then
+ write(6,*) '!!!ATTENTION i =', i, ' j = ', j, ' row = ', righe, &
+ 'col = ', colonne, 'Out of region'
+ go to 10001
+endif
+reticolo2d(i,j)=LIVELLO
+iii=iii+1
+if(iii==1) then
+ go to 1111
+endif
+
+
+!****************************************************************************************************************************************
+
+
+
+!****************************************************************************************************************************************
+!OGNI VOLTA
+10000 if(vettore_direzione>-0.414 .and. vettore_direzione<0 .and. dir=='+') then
+ go to 30001
+else if (vettore_direzione<0.414 .and. vettore_direzione>=0 .and. dir=='-') then
+ go to 30001
+else if (vettore_direzione<-0.414 .and. vettore_direzione>=-2.414 .and. dir=='+') then
+ go to 40001
+else if (abs(vettore_direzione)>2.414 .and. dir=='+') then
+ go to 50001
+else if (vettore_direzione<=2.414 .and. vettore_direzione>0.414 .and. dir=='+') then
+ go to 60001
+else if (vettore_direzione<=0.414 .and. vettore_direzione>=0 .and. dir=='+') then
+ go to 30002
+else if (vettore_direzione>-0.414 .and. vettore_direzione<0 .and. dir=='-') then
+ go to 30002
+else if (vettore_direzione<=-0.414 .and. vettore_direzione>-2.414 .and. dir=='-') then
+ go to 40002
+else if (abs(vettore_direzione)>2.414 .and. dir=='-') then
+ go to 50002
+else if (vettore_direzione<=2.414 .and. vettore_direzione>0.414 .and. dir=='-') then
+ go to 60002
+else
+ write(6,*) 'ATTENTION → problem in subroutine muoviti'
+endif
+!****************************************************************************************************************************************
+
+
+
+
+!*****************************************************************************************************************************************
+30001 vettore_direzione=4
+if(reticolo2d(i,j-1)/=LIVELLO)then
+ quota_min=quota_terreno(i,j-1)
+ vettore_direzione=0
+ dir='-'
+else
+ quota_min=LIVELLO
+endif
+if(quota_terreno(i-1,j-1)<quota_min .and. reticolo2d(i-1,j-1)/=LIVELLO)then
+ quota_min=quota_terreno(i-1,j-1)
+ vettore_direzione=-1
+ dir='+'
+endif
+if(quota_terreno(i+1,j-1)<quota_min .and. reticolo2d(i+1,j-1)/=LIVELLO)then
+ quota_min=quota_terreno(i+1,j-1)
+ vettore_direzione=1
+ dir='-'
+endif
+if(quota_terreno(i-1,j)<quota_min .and. reticolo2d(i-1,j)/=LIVELLO)then
+ quota_min=quota_terreno(i-1,j)
+ vettore_direzione=3
+ dir='+'
+endif
+if(quota_terreno(i+1,j)<quota_min .and. reticolo2d(i+1,j)/=LIVELLO)then
+ quota_min=quota_terreno(i+1,j)
+ vettore_direzione=3
+ dir='-'
+endif
+go to 70000
+
+!****************************************************************************************************************************************
+30002 vettore_direzione=4
+if(reticolo2d(i,j+1)/=LIVELLO)then
+ quota_min=quota_terreno(i,j+1)
+ vettore_direzione=0
+ dir='+'
+else
+ quota_min=LIVELLO
+endif
+if(quota_terreno(i-1,j+1)<quota_min .and. reticolo2d(i-1,j+1)/=LIVELLO)then
+ quota_min=quota_terreno(i-1,j+1)
+ vettore_direzione=1
+ dir='+'
+endif
+if(quota_terreno(i+1,j+1)<quota_min .and. reticolo2d(i+1,j+1)/=LIVELLO)then
+ quota_min=quota_terreno(i+1,j+1)
+ vettore_direzione=-1
+ dir='-'
+endif
+if(quota_terreno(i-1,j)<quota_min .and. reticolo2d(i-1,j)/=LIVELLO)then
+ quota_min=quota_terreno(i-1,j)
+ vettore_direzione=3
+ dir='+'
+endif
+if(quota_terreno(i+1,j)<quota_min .and. reticolo2d(i+1,j)/=LIVELLO)then
+ quota_min=quota_terreno(i+1,j)
+ vettore_direzione=3
+ dir='-'
+endif
+
+go to 70000
+
+
+!****************************************************************************************************************************************
+40001 vettore_direzione=4
+if(reticolo2d(i-1,j-1)/=LIVELLO)then
+ quota_min=quota_terreno(i-1,j-1)
+ vettore_direzione=-1
+ dir='+'
+else
+ quota_min=LIVELLO
+endif
+if(quota_terreno(i-1,j)<quota_min .and. reticolo2d(i-1,j)/=LIVELLO)then
+ quota_min=quota_terreno(i-1,j)
+ vettore_direzione=3
+ dir='+'
+endif
+if(quota_terreno(i,j-1)<quota_min .and. reticolo2d(i,j-1)/=LIVELLO)then
+ quota_min=quota_terreno(i,j-1)
+ vettore_direzione=0
+ dir='-'
+endif
+if(quota_terreno(i-1,j+1)<quota_min .and. reticolo2d(i-1,j+1)/=LIVELLO)then
+ quota_min=quota_terreno(i-1,j+1)
+ vettore_direzione=1
+ dir='+'
+endif
+if(quota_terreno(i+1,j-1)<quota_min .and. reticolo2d(i+1,j-1)/=LIVELLO)then
+ quota_min=quota_terreno(i+1,j-1)
+ vettore_direzione=1
+ dir='-'
+endif
+
+go to 70000
+
+!****************************************************************************************************************************************
+40002 vettore_direzione=4
+if(reticolo2d(i+1,j+1)/=LIVELLO)then
+ quota_min=quota_terreno(i+1,j+1)
+ vettore_direzione=-1
+ dir='-'
+else
+ quota_min=LIVELLO
+endif
+if(quota_terreno(i,j+1)<quota_min .and. reticolo2d(i,j+1)/=LIVELLO)then
+ quota_min=quota_terreno(i,j+1)
+ vettore_direzione=0
+ dir='+'
+endif
+if(quota_terreno(i+1,j)<quota_min .and. reticolo2d(i+1,j)/=LIVELLO)then
+ quota_min=quota_terreno(i+1,j)
+ vettore_direzione=3
+ dir='-'
+endif
+if(quota_terreno(i-1,j+1)<quota_min .and. reticolo2d(i-1,j+1)/=LIVELLO)then
+ quota_min=quota_terreno(i-1,j+1)
+ vettore_direzione=1
+ dir='+'
+endif
+if(quota_terreno(i+1,j-1)<quota_min .and. reticolo2d(i+1,j-1)/=LIVELLO)then
+ quota_min=quota_terreno(i+1,j-1)
+ vettore_direzione=1
+ dir='-'
+endif
+go to 70000
+
+!****************************************************************************************************************************************
+50001 vettore_direzione=4
+if(reticolo2d(i-1,j)/=LIVELLO)then
+ quota_min=quota_terreno(i-1,j)
+ vettore_direzione=3
+ dir='+'
+else
+ quota_min=LIVELLO
+endif
+if(quota_terreno(i-1,j+1)<quota_min .and. reticolo2d(i-1,j+1)/=LIVELLO)then
+ quota_min=quota_terreno(i-1,j+1)
+ vettore_direzione=1
+ dir='+'
+endif
+if(quota_terreno(i-1,j-1)<quota_min .and. reticolo2d(i-1,j-1)/=LIVELLO)then
+ quota_min=quota_terreno(i-1,j-1)
+ vettore_direzione=-1
+ dir='+'
+endif
+if(quota_terreno(i,j+1)<quota_min .and. reticolo2d(i,j+1)/=LIVELLO)then
+ quota_min=quota_terreno(i,j+1)
+ vettore_direzione=0
+ dir='+'
+endif
+if(quota_terreno(i,j-1)<quota_min .and. reticolo2d(i,j-1)/=LIVELLO)then
+ quota_min=quota_terreno(i,j-1)
+ vettore_direzione=0
+ dir='-'
+endif
+go to 70000
+!****************************************************************************************************************************************
+50002 vettore_direzione=4
+if(reticolo2d(i+1,j)/=LIVELLO)then
+ quota_min=quota_terreno(i+1,j)
+ vettore_direzione=3
+ dir='-'
+else
+ quota_min=LIVELLO
+endif
+if(quota_terreno(i+1,j+1)<quota_min .and. reticolo2d(i+1,j+1)/=LIVELLO)then
+ quota_min=quota_terreno(i+1,j+1)
+ vettore_direzione=-1
+ dir='-'
+endif
+if(quota_terreno(i+1,j-1)<quota_min .and. reticolo2d(i+1,j-1)/=LIVELLO)then
+ quota_min=quota_terreno(i+1,j-1)
+ vettore_direzione=1
+ dir='-'
+endif
+if(quota_terreno(i,j+1)<quota_min .and. reticolo2d(i,j+1)/=LIVELLO)then
+ quota_min=quota_terreno(i,j+1)
+ vettore_direzione=0
+ dir='+'
+endif
+if(quota_terreno(i,j-1)<quota_min .and. reticolo2d(i,j-1)/=LIVELLO)then
+ quota_min=quota_terreno(i,j-1)
+ vettore_direzione=0
+ dir='-'
+endif
+go to 70000
+!****************************************************************************************************************************************
+60001 vettore_direzione=4
+if(reticolo2d(i-1,j+1)/=LIVELLO)then
+ quota_min=quota_terreno(i-1,j+1)
+ vettore_direzione=1
+ dir='+'
+else
+ quota_min=LIVELLO
+endif
+if(quota_terreno(i,j+1)<quota_min .and. reticolo2d(i,j+1)/=LIVELLO)then
+ quota_min=quota_terreno(i,j+1)
+ vettore_direzione=0
+ dir='+'
+endif
+if(quota_terreno(i-1,j)<quota_min .and. reticolo2d(i-1,j)/=LIVELLO)then
+ quota_min=quota_terreno(i-1,j)
+ vettore_direzione=3
+ dir='+'
+endif
+if(quota_terreno(i+1,j+1)<quota_min .and. reticolo2d(i+1,j+1)/=LIVELLO)then
+ quota_min=quota_terreno(i+1,j+1)
+ vettore_direzione=-1
+ dir='-'
+endif
+if(quota_terreno(i-1,j-1)<quota_min .and. reticolo2d(i-1,j-1)/=LIVELLO)then
+ quota_min=quota_terreno(i-1,j-1)
+ vettore_direzione=-1
+ dir='+'
+endif
+go to 70000
+!****************************************************************************************************************************************
+60002 vettore_direzione=4
+if(reticolo2d(i+1,j-1)/=LIVELLO)then
+ quota_min=quota_terreno(i+1,j-1)
+ vettore_direzione=1
+ dir='-'
+else
+ quota_min=LIVELLO
+endif
+if(quota_terreno(i+1,j)<quota_min .and. reticolo2d(i+1,j)/=LIVELLO)then
+ quota_min=quota_terreno(i+1,j)
+ vettore_direzione=3
+ dir='-'
+endif
+if(quota_terreno(i,j-1)<quota_min .and. reticolo2d(i,j-1)/=LIVELLO)then
+ quota_min=quota_terreno(i,j-1)
+ vettore_direzione=0
+ dir='-'
+endif
+if(quota_terreno(i+1,j+1)<quota_min .and. reticolo2d(i+1,j+1)/=LIVELLO)then
+ quota_min=quota_terreno(i+1,j+1)
+ vettore_direzione=-1
+ dir='-'
+endif
+if(quota_terreno(i-1,j-1)<quota_min .and. reticolo2d(i-1,j-1)/=LIVELLO)then
+ quota_min=quota_terreno(i-1,j-1)
+ vettore_direzione=-1
+ dir='+'
+endif
+
+
+
+70000 if (vettore_direzione==1) then
+ if(dir=='+')then
+ i=i-1
+ j=j+1
+ else if (dir=='-')then
+ i=i+1
+ j=j-1
+ endif
+else if (vettore_direzione==3) then
+ if(dir=='+')then
+ i=i-1
+ j=j
+ else if (dir=='-')then
+ i=i+1
+ j=j
+ endif
+else if (vettore_direzione==0) then
+ if(dir=='+')then
+ i=i
+ j=j+1
+ else if (dir=='-')then
+ i=i
+ j=j-1
+ endif
+else if (vettore_direzione==-1) then
+ if(dir=='+')then
+ i=i-1
+ j=j-1
+ else if (dir=='-')then
+ i=i+1
+ j=j+1
+ endif
+endif
+
+
+
+
+
+!****************************************************************************************************************************************
+!CONTROLLI DI ROUTINE ----> se li supero torno a 100000
+! ----> se non li supero ho finito di muovermi
+ if (quota_terreno(i,j)>=LIVELLO) then
+ reticolo2d(i,j)=0
+ go to 10001
+else if (i>=righe .or. j>=colonne)then
+ reticolo2d(i,j)=0
+ go to 10001
+else if ( bordi_alveo(i,j)==1)then
+ go to 10001
+else if (reticolo2d(i,j)>LIVELLO)then
+ go to 10001
+else if (quota_terreno(i,j)==0) then
+ reticolo2d(i,j)=0
+ go to 10001
+else if (i==1 .or. j==1 .or. i==righe .or. j==colonne)then
+ reticolo2d(i,j)=LIVELLO
+ go to 10001
+else
+!else if (inondazioneT(i,j)==0)then
+ iiii=i
+ jjjj=j
+ !write(6,*)'controllo con il 2x2', vettore_direzione, dir
+ go to 20000
+!else
+! reticolo2d(i,j)=LIVELLO
+! go to 10000
+endif
+
+
+!ripercorro l'ultimo tratto con risoluzione del dtm 2x2 per controllare che l'acqua non abbia superato un argine
+20000 passo=res2
+if(vettore_direzione==1 .and. dir=='+') then
+ i=i+1
+ j=j-1
+ E_t=ovest+res*j-res/2
+ N_t=nord-(res*i-res/2)
+ !write(6,*)est, E_t, ovest, sud, N_t, nord
+ xxxx=0
+ 11 xxxx=xxxx+passo
+ E_t=E_t+passo/sqrt(2.)
+ N_t=N_t+passo/sqrt(2.)
+ if(xxxx>(sqrt(2.)*res)) then
+ i=i-1
+ j=j+1
+ reticolo2d(iiii,jjjj)=LIVELLO
+ go to 10000
+ endif
+ i_temp=1+int((nord-N_t)/res2)
+ j_temp=1+int((E_t-ovest)/res2)
+ !write(6,*) i_temp, righe2, J_temp,colonne2, inondazioneT_2(i_temp,j_temp)
+ call trova_elemento(i_temp,j_temp)
+ if(output>livello)then
+ ! write(6,*) '0'
+ reticolo2d(iiii,jjjj)=0
+ go to 10001
+ else
+ !write(6,*) '1'
+ go to 11
+ endif
+else if (vettore_direzione==0 .and. dir=='+' ) then
+ i=i
+ j=j-1
+ E_t=ovest+res*j-res/2
+ N_t=nord-(res*i-res/2)
+ ! write(6,*)est, E_t, ovest, sud, N_t, nord
+ xxxx=0
+ 15 xxxx=xxxx+passo
+ E_t=E_t+passo
+ if(xxxx>(res)) then
+ i=i
+ j=j+1
+ reticolo2d(iiii,jjjj)=LIVELLO
+ go to 10000
+ endif
+ i_temp=1+int((nord-N_t)/res2)
+ j_temp=1+int((E_t-ovest)/res2)
+ call trova_elemento(i_temp,j_temp)
+ if(output>LIVELLO)then
+ !write(6,*) '0'
+ reticolo2d(iiii,jjjj)=0
+ go to 10001
+ else
+ !write(6,*) '1'
+ go to 15
+ endif
+else if (vettore_direzione==-1 .and. dir=='-') then
+ i=i-1
+ j=j-1
+ E_t=ovest+res*j-res/2
+ N_t=nord-(res*i-res/2)
+ xxxx=0
+ 12 xxxx=xxxx+passo
+ E_t=E_t+passo/sqrt(2.)
+ N_t=N_t-passo/sqrt(2.)
+ if(xxxx>(sqrt(2.)*res)) then
+ j=j+1
+ i=i+1
+ reticolo2d(iiii,jjjj)=LIVELLO
+ go to 10000
+ endif
+ i_temp=1+int((nord-N_t)/res2)
+ j_temp=1+int((E_t-ovest)/res2)
+ call trova_elemento(i_temp,j_temp)
+ if(output>LIVELLO)then
+ reticolo2d(iiii,jjjj)=0
+ go to 10001
+ else
+ go to 12
+ endif
+else if (vettore_direzione==3 .and. dir=='-') then
+ i=i-1
+ j=j
+ E_t=ovest+res*j-res/2
+ N_t=nord-(res*i-res/2)
+ xxxx=0
+ 17 xxxx=xxxx+passo
+ N_t=N_t-passo
+ if(xxxx>res) then
+ i=i+1
+ j=j
+ reticolo2d(iiii,jjjj)=LIVELLO
+ go to 10000
+ endif
+ i_temp=1+int((nord-N_t)/res2)
+ j_temp=1+int((E_t-ovest)/res2)
+ call trova_elemento(i_temp,j_temp)
+ if(output>LIVELLO)then
+ reticolo2d(iiii,jjjj)=0
+ go to 10001
+ else
+ go to 17
+ endif
+else if (vettore_direzione==1 .and. dir=='-') then
+ i=i-1
+ j=j+1
+ E_t=ovest+res*j-res/2
+ N_t=nord-(res*i-res/2)
+ xxxx=0
+ 13 xxxx=xxxx+passo
+ E_t=E_t-passo/sqrt(2.)
+ N_t=N_t-passo/sqrt(2.)
+ if(xxxx>(sqrt(2.)*res)) then
+ i=i+1
+ j=j-1
+ reticolo2d(iiii,jjjj)=LIVELLO
+ go to 10000
+ endif
+ i_temp=1+int((nord-N_t)/res2)
+ j_temp=1+int((E_t-ovest)/res2)
+ call trova_elemento(i_temp,j_temp)
+ if(output>LIVELLO)then
+ reticolo2d(iiii,jjjj)=0
+ go to 10001
+ else
+ go to 13
+ endif
+else if (vettore_direzione==0 .and. dir=='-') then
+ i=i
+ j=j+1
+ E_t=ovest+res*j-res/2
+ N_t=nord-(res*i-res/2)
+ xxxx=0
+ 16 xxxx=xxxx+passo
+ E_t=E_t-passo
+ if(xxxx>(res)) then
+ i=i
+ j=j-1
+ reticolo2d(iiii,jjjj)=LIVELLO
+ go to 10000
+ endif
+ i_temp=1+int((nord-N_t)/res2)
+ j_temp=1+int((E_t-ovest)/res2)
+ call trova_elemento(i_temp,j_temp)
+ if(output>LIVELLO)then
+ reticolo2d(iiii,jjjj)=0
+ go to 10001
+ else
+ go to 16
+ endif
+else if (vettore_direzione==-1 .and. dir=='+') then
+ i=i+1
+ j=j+1
+ E_t=ovest+res*j-res/2
+ N_t=nord-(res*i-res/2)
+ xxxx=0
+ 14 xxxx=xxxx+passo
+ E_t=E_t-passo/sqrt(2.)
+ N_t=N_t+passo/sqrt(2.)
+ if(xxxx>(sqrt(2.)*res)) then
+ i=i-1
+ j=j-1
+ reticolo2d(iiii,jjjj)=LIVELLO
+ go to 10000
+ endif
+ i_temp=1+int((nord-N_t)/res2)
+ j_temp=1+int((E_t-ovest)/res2)
+ call trova_elemento(i_temp,j_temp)
+ if(output>LIVELLO)then
+ reticolo2d(iiii,jjjj)=0
+ go to 10001
+ else
+ go to 14
+ endif
+else if (vettore_direzione==3 .and. dir=='+') then
+ i=i+1
+ j=j
+ E_t=ovest+res*j-res/2
+ N_t=nord-(res*i-res/2)
+ xxxx=0
+ 18 xxxx=xxxx+passo
+ N_t=N_t+passo
+ if(xxxx>res) then
+ i=i-1
+ j=j
+ reticolo2d(iiii,jjjj)=LIVELLO
+ go to 10000
+ endif
+ i_temp=1+int((nord-N_t)/res2)
+ j_temp=1+int((E_t-ovest)/res2)
+ call trova_elemento(i_temp,j_temp)
+ if(output>LIVELLO)then
+ reticolo2d(iiii,jjjj)=0
+ go to 10001
+ else
+ go to 18
+ endif
+endif
+
+
+10001 temp=0
+
+end subroutine muoviti
+
+
+
+
+
+
+
+
+
+!************************************************************************************************
+!Subroutine allocazione matrice sparsa
+!************************************************************************************************
+
+subroutine matrice_sparsa (nomefile,righe,colonne)
+!uso l'algoritmo CSR(Compressed )
+use dd
+implicit none
+character(40)::nomefile
+integer::i,j, cont, righe,colonne, nn, cont2
+real::a
+real,allocatable,dimension(:):: temp
+
+
+allocate(temp(colonne))
+allocate (iind(righe+1))
+open (unit=70, file=nomefile)
+cont=0
+
+do i=1,righe
+ read(70,*)(temp(j), j=1,colonne)
+ !write(6,*) temp
+ do j=1,colonne
+ if (temp(j)/=0.) then
+ cont=cont+1
+ endif
+ enddo
+enddo
+rewind(70)
+
+n_no_null=cont
+write(6,*)'There are', n_no_null, 'no-null cells in the Digital Terrain Model'
+
+allocate (val(n_no_null))
+allocate (jind(n_no_null))
+
+
+cont=0
+cont2=0
+do i=1,righe
+ cont2=cont+1
+ iind(i)=cont2
+ read (70,*) (temp(j),j=1,colonne)
+ do j=1,colonne
+ if (temp(j)/=0.) then
+ cont=cont+1
+ val(cont)=temp(j)
+ jind(cont)=j
+ endif
+ enddo
+enddo
+iind(righe+1)=n_no_null+1
+
+!write(6,*) "val=", val
+!write(6,*) "iind=", iind
+!write(6,*) "jind=", jind
+
+end subroutine matrice_sparsa
+
+
+
+
+
+!************************************************************************************************
+!Subroutine trova punto in matrice sparsa
+!************************************************************************************************
+
+subroutine trova_elemento(iriga,jcolonna)
+use dd
+integer::i,iriga,jcolonna, iniz, fine
+
+iniz=iind(iriga)
+fine=iind(iriga+1)-1
+
+do i=iniz,fine
+ if (jind(i)==jcolonna)then
+ output=val(i)
+ go to 144
+ else
+ output=0.
+ endif
+enddo
+
+144 i=0
+
+end subroutine trova_elemento
+
+
+
Added: grass-addons/raster/r.inund.fluv/README.txt
===================================================================
--- grass-addons/raster/r.inund.fluv/README.txt (rev 0)
+++ grass-addons/raster/r.inund.fluv/README.txt 2008-05-02 12:24:52 UTC (rev 31204)
@@ -0,0 +1,41 @@
+This is a script written in bash shell language, configured for grass 6.3 release.
+Users of previous releases can use the r.inund.fluv_62 script (see the notes below)
+
+r.inund.fluv use the following fortran codes:
+- find_main_channel.f90
+- clean_inundation.f90
+- 2d_path.f90
+- correction_from_path.f90
+
+***************************REQUEST ********************************
+The installation of this code requests a fortran compiler installed on your computer.
+The installation script (install.sh) is written for gnu fortran compiler with optimization.
+(http://www.gnu.org/software/gcc/fortran/)
+If you have an other compiler you may change the install script (install.sh)
+
+e.g.
+gfortran -O1 -o ***.exe ***.f90
+
+where
+*) `gfortran´ is the name of your compiler
+*) `-O1´ is the optimization option
+*) `-o ***.exe´ is the name of output executable (DO NOT CHANGE!)
+*) `***.f90´ is the name of code (DO NOT CHANGE!)
+
+For any other informations contact the autors of the script
+(see description.html)
+******************************************************************
+
+**************** Users of GRASS 6.3 ******************************
+1) Run the install.sh script
+sh install.sh
+
+
+**************** Users of GRASS 6.2 ******************************
+1) Change the name of the script for GRASS 6.3.x releases
+mv r.inund.fluv r.inund.fluv_63
+2) Change the name of script for grass 6.2.x releases
+mv r.inund.fluv_62 r.inund.fluv
+3) Run the install.sh script
+sh install.sh
+
Added: grass-addons/raster/r.inund.fluv/clean_inundation.f90
===================================================================
--- grass-addons/raster/r.inund.fluv/clean_inundation.f90 (rev 0)
+++ grass-addons/raster/r.inund.fluv/clean_inundation.f90 2008-05-02 12:24:52 UTC (rev 31204)
@@ -0,0 +1,376 @@
+module dd
+integer,parameter::ferma_bisezione=15
+integer::n_no_null
+real, allocatable, dimension(:)::val, iind, jind
+real::output
+end module dd
+
+
+
+!*********************************************************************************************************************************
+!il programma si muove su ogni pixel della carta raster di inondazione che è già stata pulita dai "laghi",
+!dove il pixel vale 0 (NON inondato) lasci stare tutto com'è
+!dove il pixel vale 1 (inondato)--> cerca il punto del profilo più vicino,
+! calcola la retta che unisce il pixel al profilo
+! si muove su questa retta e controlla se i pixel attraversati sono asciutti o bagnati
+! se asciutti--> il pixel diventa anch'esso asciutto
+! se bagnati il pixel rimane bagnato
+!*********************************************************************************************************************************
+
+program pulizia
+use dd
+implicit none
+real::est, ovest, nord, sud, nord1, sud1, est1, ovest1
+integer::i, j, ii, righe, colonne, iend, nrighe, int_x, int_y, row, col
+integer::pixel_i, pixel_j, temp_x, temp_y, righe1, colonne1
+integer, allocatable, dimension(:,:):: inond_nuova, inond
+real(kind=8)::passo, ver, coord_est, coord_nord, x, y, dist, dist_temp, E1, N1,res, res1
+real, allocatable,dimension(:)::E, N, q
+character(len=20)::xx, xx1
+real, allocatable, dimension(:,:)::quota
+character(40)::nomefile1, nomefile2, nomefile3, nomefile4
+character(40)::nomefile5, nomefile6, nomefile7
+open(unit=40, status='old', file='nomefile.txt')
+read(40,*)nomefile1, nomefile2, nomefile3, nomefile4, nomefile5, nomefile6, nomefile7
+close(40)
+
+
+
+!apertura file
+!********************************************************************************************************************************
+!N.B. controlla il nome file
+!********************************************************************************************************************************
+!50=region20
+!51=quote_punti_adiacenti (superficie idrica ottenuta con interpolazone lineare)
+!52=inondazione_step2
+!53=file con il profilo
+!54=region2
+!dtm2x2= nomefile7 (vd dopo)
+!60=inondazione_nuova (step3)
+
+open (unit=50, status='old', file=nomefile5)
+open (unit=51, status='old', file=nomefile4)
+open (unit=52, status='old', file=nomefile1)
+open (unit=53, status='old', file=nomefile2)
+open (unit=54, status='old', file=nomefile6)
+!open (unit=55, status='old', file='grass_script/dtm2x2')
+open (unit=60, file=nomefile3)
+
+
+
+!guardo quant'è lungo il file che contiene il profilo
+do i=1,10000000
+ read(53,*,iostat=iend)xx
+ if(iend /= 0)exit
+enddo
+if(iend > 0)then
+ write(6,*)'Error reading file of water surface pofile:',iend
+ stop
+endif
+nrighe=i-1
+write(6,*) 'The file of water surface profile has ',nrighe,' rows'
+rewind(53)
+
+allocate(E(nrighe))
+allocate(N(nrighe))
+allocate(q(nrighe))
+
+do i=1,nrighe
+ read(53,*) E(i), N(i), q(i)
+enddo
+
+!leggo le dimensioni della regione meno dettagliata
+do i=1,13
+ read(50,*)xx
+ if (xx=='north:') then
+ backspace(50)
+ read(50,*)xx, nord
+ ! write(60,'(a6,i)') trim(xx), nord
+ write(60,'(a6,f10.2)') trim(xx), nord
+ !write(60,*) trim(xx), nord
+ else if (xx=='south:') then
+ backspace(50)
+ read(50,*)xx, sud
+ write(60,'(a6,f10.2)')trim(xx), sud
+ !write(60,*) trim(xx), sud
+ else if (xx=='east:') then
+ backspace(50)
+ read(50,*)xx, est
+ write(60,'(a5,f10.2)') trim(xx), est
+ !write(60,*) trim(xx), est
+ else if (xx=='west:') then
+ backspace(50)
+ read(50,*)xx, ovest
+ write(60,'(a5,f10.2)')trim(xx), ovest
+ !write(60,*) trim(xx), ovest
+ else if (xx=='rows:') then
+ backspace(50)
+ read(50,*)xx, righe
+ write(60,'(a5,i5)')trim(xx), righe
+ !write(60,*) trim(xx), righe
+ else if (xx=='cols:') then
+ backspace(50)
+ read(50,*)xx, colonne
+ write(60,'(a5,i5)')trim(xx), colonne
+ !write(60,*) trim(xx), colonne
+ else if (xx=='nsres:') then
+ backspace(50)
+ read(50,*)xx, res
+ endif
+enddo
+write(6,*) nord, sud, est, ovest, righe, colonne
+close(50)
+
+
+
+
+!alloco le matrici del file
+allocate(quota(righe,colonne))
+allocate(inond_nuova(righe,colonne))
+allocate(inond(righe,colonne))
+
+read(51,*) ((quota(i,j), j=1,colonne), i=1,righe)
+write(6,*) 'Read the raster map with water surface'
+read(52,*) ((inond(i,j), j=1,colonne), i=1,righe)
+write(6,*) 'Read the raster inondation map of step 1'
+close(51)
+close(52)
+
+!leggo le dimensioni della regione con risoluzione dettagliata
+do i=1,13
+ read(54,*)xx
+ if (xx=='north:') then
+ backspace(54)
+ read(54,*)xx, nord1
+ else if (xx=='south:') then
+ backspace(54)
+ read(54,*)xx, sud1
+ else if (xx=='east:') then
+ backspace(54)
+ read(54,*)xx, est1
+ else if (xx=='west:') then
+ backspace(54)
+ read(54,*)xx, ovest1
+ else if (xx=='rows:') then
+ backspace(54)
+ read(54,*)xx, righe1
+ else if (xx=='cols:') then
+ backspace(54)
+ read(54,*)xx, colonne1
+ else if (xx=='nsres:') then
+ backspace(54)
+ read(54,*)xx, res1
+ endif
+enddo
+write(6,*) nord1, sud1, est1, ovest1, righe1, colonne1, res1
+
+
+!open(unit=70,file='control.txt')
+!write(70,*)res1
+!close(70)
+!allocate(quota_terreno(righe1,colonne1))
+
+!leggo il DTM
+!read(55,*) ((quota_terreno(i,j), j=1,colonne1), i=1,righe1)
+!write(6,*) 'letto il DTM'
+!rewind(55)
+!close(55)
+
+call matrice_sparsa (nomefile7,righe1,colonne1)
+write(6,*) 'Read the Digital Terrain Model with Compressed Sparse Row Algorithm'
+
+!scorro tutti i pixel della matrice inondazione
+!se trovo 0 (area NON inondata) ricopio il valore nella nuova matrice
+!se trovo 1 eseguo il controllo e vedo se lasciare 1 o 0
+
+j=0
+i=0
+1000 j=j+1
+1001 i=i+1
+!write(6,*)i,j
+if(i==righe+1) then
+ if(j==colonne)then
+ go to 3000
+ endif
+ i=0
+ go to 1000
+endif
+
+if(inond(i,j)==0) then
+ inond_nuova(i,j)=inond(i,j)
+ go to 1001
+else if(inond(i,j)==1) then
+ coord_est=ovest+res*j-res/2
+ coord_nord=nord-(res*i-res/2)
+ !write(6,*) coord_est,coord_nord,i,j
+ !cerco il punto del profilo più vicino al mio pixel
+ dist=100000
+ do ii=1,nrighe
+ if (q(ii)==quota(i,j)) then
+ dist_temp=sqrt((coord_est-E(ii))**2+(coord_nord-N(ii))**2)
+ if(dist_temp<=dist) then
+ E1=E(ii)
+ N1=N(ii)
+ dist=dist_temp
+ endif
+ endif
+ enddo
+ !controllo
+ if(E1>est .or. E1<ovest .or. N1>nord .or. N1<sud) then
+ write(6, *) 'A point of water surface profile is out of region Est=',E1, 'Nord=',N1
+ if(E1>est)then
+ E1=est
+ dist=sqrt((coord_est-E1)**2+(coord_nord-N1)**2)
+ else if (E1<ovest) then
+ E1=ovest
+ dist=sqrt((coord_est-E1)**2+(coord_nord-N1)**2)
+ else if (N1>nord) then
+ N1=nord
+ dist=sqrt((coord_est-E1)**2+(coord_nord-N1)**2)
+ else if (N1<sud) then
+ N1=sud
+ dist=sqrt((coord_est-E1)**2+(coord_nord-N1)**2)
+ endif
+ endif
+ rewind(53)
+
+ !mi muovo sulla retta che congiunge l'alveo al pixel con un passo pari alla risoluzione più dettagliata
+ passo=res1
+ ver=0.
+ x=coord_est
+ y=coord_nord
+ 7000 ver=ver+passo
+ x=x+passo*(E1-coord_est)/dist
+ y=y+passo*(N1-coord_nord)/dist
+ !ora devo dire in quale pixel sono
+ pixel_i=1+int((nord1-y)/res1)
+ pixel_j=1+int((x-ovest1)/res1)
+ if (ver>dist) then
+ inond_nuova(i,j)=inond(i,j)
+ go to 1001
+ endif
+ call trova_elemento(pixel_i,pixel_j)
+ if (output>quota(i,j))then
+ inond_nuova(i,j)=0
+ go to 1001
+ else
+ go to 7000
+ endif
+endif
+
+
+
+!scrivo la nuova matrice
+3000 write(6,*) 'Writing output raster Map'
+222 format(5000(i1,1x))
+write(60,222) ((inond_nuova(i,j), j=1,colonne),i=1,righe)
+!write(60,'(i1)') ((inond_nuova(i,j), j=1,colonne),i=1,righe)
+!write(61,'(b)') ((inond_nuova(i,j), j=1,colonne),i=1,righe)
+
+
+deallocate(E)
+deallocate(N)
+deallocate(q)
+deallocate(quota)
+deallocate(inond_nuova)
+deallocate(inond)
+
+
+write(6,*)"END of Fortran code → OK"
+end program pulizia
+
+
+
+
+
+
+
+
+
+!************************************************************************************************
+!Subroutine allocazione matrice sparsa
+!************************************************************************************************
+
+subroutine matrice_sparsa (nomefile,righe,colonne)
+!uso l'algoritmo CSR(Compressed )
+use dd
+implicit none
+character(40)::nomefile
+integer::i,j, cont, righe,colonne, nn, cont2
+real::a
+real,allocatable,dimension(:):: temp
+
+
+allocate(temp(colonne))
+allocate (iind(righe+1))
+open (unit=70, file=nomefile)
+cont=0
+
+do i=1,righe
+ read(70,*)(temp(j), j=1,colonne)
+ !write(6,*) temp
+ do j=1,colonne
+ if (temp(j)/=0.) then
+ cont=cont+1
+ endif
+ enddo
+enddo
+rewind(70)
+
+n_no_null=cont
+
+write(6,*)'There are', n_no_null, 'no null cells in Digital Terrain Model'
+allocate (val(n_no_null))
+allocate (jind(n_no_null))
+
+
+cont=0
+cont2=0
+do i=1,righe
+ cont2=cont+1
+ iind(i)=cont2
+ read (70,*) (temp(j),j=1,colonne)
+ do j=1,colonne
+ if (temp(j)/=0.) then
+ cont=cont+1
+ val(cont)=temp(j)
+ jind(cont)=j
+ endif
+ enddo
+enddo
+iind(righe+1)=n_no_null+1
+
+!write(6,*) "val=", val
+!write(6,*) "iind=", iind
+!write(6,*) "jind=", jind
+
+end subroutine matrice_sparsa
+
+
+
+
+
+!************************************************************************************************
+!Subroutine trova punto in matrice sparsa
+!************************************************************************************************
+
+subroutine trova_elemento(iriga,jcolonna)
+use dd
+integer::i,iriga,jcolonna, iniz, fine
+
+iniz=iind(iriga)
+fine=iind(iriga+1)-1
+
+do i=iniz,fine
+ if (jind(i)==jcolonna)then
+ output=val(i)
+ go to 144
+ else
+ output=0.
+ endif
+enddo
+
+144 i=0
+
+end subroutine trova_elemento
+
Added: grass-addons/raster/r.inund.fluv/correction_from_path.f90
===================================================================
--- grass-addons/raster/r.inund.fluv/correction_from_path.f90 (rev 0)
+++ grass-addons/raster/r.inund.fluv/correction_from_path.f90 2008-05-02 12:24:52 UTC (rev 31204)
@@ -0,0 +1,615 @@
+
+module dd
+integer,parameter::ferma_bisezione=15
+integer::n_no_null
+real, allocatable, dimension(:)::val, iind, jind
+real::output
+end module dd
+
+
+!*******************************************************************************************************************************************************
+!il programma controlla quelle aree che il precedente programma fortran aveva pulito tenendo in conto la presenza di argini
+!l'acqua infatti può arrivare al di là degli argini anche percorrendo percorsi NON ortogonali all'alveo come si è ivece ipotizzato nei primi 2 step
+!in particolare l'acqua può seguire alcuni percorsi secondo la massima pendenza---> a tale scopo si è precedentemente creato un reticolo
+!il programma per ogni pixel == 1 della matrice diff_inond_2 cerca il punto del reticolo2d più vicino
+!quindi percorre il tratto tra il pixel e il reticolo
+!se incontra DTM più alto o punto asciutto---> pixel==0
+!in caso contrario ----> pixel=1
+! a tale scopo utilizza in input: - il file contenente i limiti della regione
+! - il file contenente la matrice diff_inond
+! - il DTM
+! - i file contenente i punti del reticolo
+!fornisce in output per il profilo considerato una matrice con le aree da inondare nuovamente che poi andrà sommata alla carta di inondazione
+!*******************************************************************************************************************************************************
+
+
+program correzione_punti
+use dd
+implicit none
+
+real:: nord, sud, ovest, est
+integer,parameter::punti=15
+integer::i,j, ii, iii, f, iend, righe, colonne, pixel_temp_i, pixel_temp_j, npunti, caso, righe2,colonne2
+integer, allocatable, dimension(:,:):: diff_inond, pixel_inondati
+!real, allocatable, dimension(:,:):: punti3d
+real(kind=8)::x, temp, dist_min, dist_temp, E_pixel, N_pixel, E_t, N_t, res2, res
+real(kind=8),dimension(punti)::dist, EE, NN,qq
+real, allocatable, dimension(:):: E, N, quota
+character(len=40)::xx, tab, tab1
+character(40)::nomefile1,nomefile2,nomefile3,nomefile4,nomefile5,nomefile6
+open(unit=40, status='old', file='nomefile.txt')
+read(40,*)nomefile1,nomefile2,nomefile3,nomefile4,nomefile5,nomefile6
+close(40)
+
+
+!50 region10.txt
+!51 diff_inond
+! dtm2x2 nomefile5 (vd dopo)
+!53 punti3d_reticolo
+!54 region2.txt
+!60 correzione
+
+open(unit=50, status='old', file=nomefile4)
+open(unit=51, status='old', file=nomefile1)
+!open(unit=52, status='old', file='grass_script/dtm2x2')
+open(unit=53, status='old', file=nomefile2)
+open(unit=54, status='old', file=nomefile6)
+open(unit=60, file=nomefile3)
+
+
+!leggo le dimensioni della regione utili per allocare le matrici in cui memorizzare carta di input e di output
+do i=1,13
+ read(50,*)xx
+ if (xx=='north:') then
+ backspace(50)
+ read(50,*)xx, nord
+ write(60,'(a6,f10.2)') trim(xx), nord
+ !write(60,'(a6,i)') trim(xx), nord
+ !write(60,*) trim(xx), nord
+ else if (xx=='nsres:') then
+ backspace(50)
+ read(50,*)xx, res
+ else if (xx=='south:') then
+ backspace(50)
+ read(50,*)xx, sud
+ write(60,'(a6,f10.2)')trim(xx), sud
+ !write(60,*) trim(xx), sud
+ else if (xx=='east:') then
+ backspace(50)
+ read(50,*)xx, est
+ write(60,'(a5,f10.2)') trim(xx), est
+ !write(60,*) trim(xx), est
+ else if (xx=='west:') then
+ backspace(50)
+ read(50,*)xx, ovest
+ write(60,'(a5,f10.2)')trim(xx), ovest
+ !write(60,*) trim(xx), ovest
+ else if (xx=='rows:') then
+ backspace(50)
+ read(50,*)xx, righe
+ write(60,'(a5,i5)')trim(xx), righe
+ !write(60,*) trim(xx), righe
+ else if (xx=='cols:') then
+ backspace(50)
+ read(50,*)xx, colonne
+ write(60,'(a5,i5)')trim(xx), colonne
+ !write(60,*) trim(xx), colonne
+ endif
+enddo
+write(6,*) nord, sud, est, ovest, righe, colonne,res
+rewind(50)
+close(50)
+
+
+allocate(diff_inond(righe,colonne))
+allocate(pixel_inondati(righe,colonne))
+!allocate(punti3d(righe,colonne))
+
+!leggo la matrice diff_inond
+read(51,*) ((diff_inond(i,j), j=1,colonne), i=1,righe)
+write(6,*) 'Read raster map of areas to control'
+rewind(51)
+close(51)
+
+do i=1,13
+ read(54,*)xx
+ if (xx=='nsres:') then
+ backspace(54)
+ read(54,*)xx, res2
+ res2=int(res2)
+ else if (xx=='rows:') then
+ backspace(54)
+ read(54,*)xx, righe2
+ else if (xx=='cols:') then
+ backspace(54)
+ read(54,*)xx, colonne2
+ endif
+enddo
+write(6,*) righe2, colonne2,res2
+rewind(54)
+close(54)
+!allocate(quota_terreno(righe2,colonne2))
+
+!leggo il DTM
+!read(52,*) ((quota_terreno(i,j), j=1,colonne2), i=1,righe2)
+!write(6,*) 'letto il DTM'
+!rewind(52)
+!close(52)
+call matrice_sparsa (nomefile5,righe2,colonne2)
+write(6,*) 'Read the Digital Terrain Model with Compressed Sparse Rows algorithm'
+
+
+!scrivo 0 sulla matrice di output
+do j=1,colonne
+ do i=1,righe
+ pixel_inondati(i,j)=0
+ enddo
+enddo
+
+!guardo quant'è lungo il file che contiene i punti
+do i=1,10000000
+ read(53,*,iostat=iend)xx
+ if(iend /= 0)exit
+enddo
+if(iend > 0)then
+ write(6,*)'Error reading points of course by each river section:',iend
+ stop
+endif
+npunti=i-1
+write(6,*) 'There are',npunti,' of water course by each river section to analize'
+rewind(53)
+
+
+allocate(E(npunti))
+allocate(N(npunti))
+allocate(quota(npunti))
+
+
+call modifica_ascii(npunti)
+
+open(unit=55,status='old',file='punti_mod')
+
+do i=1,npunti
+ read(55,*) E(i), N(i), quota(i)
+ !write(6,*) i, E(i), N(i), quota(i)
+enddo
+!write(6,*) 'fatto per oggi basta!!!'
+
+
+j=0
+i=0
+1000 j=j+1
+1001 i=i+1
+!write(6,*)i,j
+if(i==righe+1) then
+ if(j==colonne)then
+ go to 1111
+ endif
+ i=0
+ go to 1000
+endif
+if(diff_inond(i,j)==1) then
+ !write(6,*)i,j
+ E_pixel=ovest+res*j-res/2
+ N_pixel=nord-(res*i-res/2)
+ dist_min=1000000000000000000.
+ !bisogna cercare i punti (il numero è 1 paramatro del programma) più vicini a detto pixel
+ !quindi si esegue un controllo sul percorso pixel punto
+ !se la quota del dtm è sempre inferiore a quella idraulica -->pixel a rischio
+ !se la quota del dtm fosse maggiore --> si passa a vedere il percorso dal punto "i+1 -esimo"
+ !conclusi i punti se il dtm è sempre + elevato ---> pixel non a rischio
+ do ii=1,punti-1
+ dist(ii)=0.
+ EE(ii)=E(ii)
+ NN(ii)=N(ii)
+ qq(ii)=quota(ii)
+ enddo
+ do ii=1,npunti
+ dist_temp=sqrt((E(ii)-E_pixel)**2.+(N(ii)-N_pixel)**2.)
+ if (dist_temp<dist_min)then
+ do iii=punti,2,-1
+ dist(iii)=dist(iii-1)
+ EE(iii)=EE(iii-1)
+ NN(iii)=NN(iii-1)
+ qq(iii)=qq(iii-1)
+ enddo
+ dist(1)=dist_temp
+ EE(1)=E(ii)
+ NN(1)=N(ii)
+ qq(1)=quota(ii)
+ dist_min=dist_temp
+ endif
+ enddo
+ !se la distanza è < della semidiagonale del pixel allora vuol dire che siamo nel pixel dove si trova il reticolo --> a rischio
+ if (dist(1)<res*sqrt(2.)/2) then
+ pixel_inondati(i,j)=1
+ go to 1001
+ endif
+ x=res2
+ caso=0
+1003 caso=caso+1
+ E_t=EE(caso)
+ N_t=NN(caso)
+ ! write(6,*)E_t, N_t,x
+ temp=0.
+ 1002 temp=temp+x
+ !write(6,*)temp, dist(caso)
+ E_t=E_t+x*(E_pixel-EE(caso))/dist(caso)
+ N_t=N_t+x*(N_pixel-NN(caso))/dist(caso)
+ pixel_temp_i=1+int((nord-N_t)/res2)
+ pixel_temp_j=1+int((E_t-ovest)/res2)
+ if(temp>dist(caso))then
+ pixel_inondati(i,j)=1
+ !write(6,*)'1', i,j
+ go to 1001
+ endif
+ ! if (pixel_temp_j==0)then
+ ! write(6,*) i, j, caso,x, temp, dist(caso), E_t, EE(caso),E_pixel, N_t, NN(caso), N_pixel
+ ! stop
+ ! endif
+ call trova_elemento(pixel_temp_i,pixel_temp_j)
+ if(qq(caso)<=output)then
+ pixel_inondati(i,j)=0
+ !write(6,*)i,j,quota_terreno(pixel_temp_i,pixel_temp_j),quota(caso),'0'
+ if(caso<punti)then
+ go to 1003
+ endif
+ go to 1001
+ else
+ go to 1002
+ endif
+else
+ go to 1001
+endif
+
+1111 write(6,*)'End of computational steps → writing output raster map'
+
+
+!scrivo la matrice dei risultati
+write(60,*) ((pixel_inondati(i,j), j=1,colonne), i=1,righe)
+
+write(6,*)'END of Fortran Code → OK'
+end program correzione_punti
+
+
+
+
+
+
+!*******************************************************************************************************************************************************
+!subroutine per modificare il file ascii
+!*******************************************************************************************************************************************************
+subroutine modifica_ascii(ferma)
+implicit none
+integer::i
+integer,intent(IN)::ferma
+character(1)::v
+character(32)::vvv
+
+write(6,*)ferma
+
+open(unit=54,file='punti_mod')
+
+i=0
+10 i=i+1
+if(i==ferma+1)then
+ write(6,*) 'Stop', ferma
+ go to 11
+endif
+read(53,*)vvv
+!write(6,*)vvv,i
+if (vvv(7:7)=='|')then
+ if (vvv(15:15)=='|') then
+ if(vvv(22:22)=='|')then
+ write(54,*)vvv(1:6), ' ', vvv(8:14), ' ', vvv(16:21)
+ go to 10
+ else if (vvv(17:17)=='|')then
+ write(54,*)vvv(1:6), ' ', vvv(8:14), ' ', vvv(16:16)
+ go to 10
+ else if (vvv(18:18)=='|')then
+ write(54,*)vvv(1:6), ' ', vvv(8:14), ' ', vvv(16:17)
+ go to 10
+ else if (vvv(19:19)=='|')then
+ write(54,*)vvv(1:6), ' ', vvv(8:14), ' ', vvv(16:18)
+ go to 10
+ else if (vvv(20:20)=='|') then
+ write(54,*)vvv(1:6), ' ', vvv(8:14), ' ', vvv(16:19)
+ go to 10
+ else if (vvv(21:21)=='|') then
+ write(54,*)vvv(1:6), ' ', vvv(8:14), ' ', vvv(16:20)
+ go to 10
+ else if (vvv(23:23)=='|') then
+ write(54,*)vvv(1:6), ' ', vvv(8:14), ' ', vvv(16:22)
+ go to 10
+ endif
+ else if(vvv(17:17)=='|') then
+ if (vvv(24:24)=='|')then
+ write(54,*)vvv(1:6), ' ', vvv(8:16), ' ', vvv(18:23)
+ go to 10
+ else if (vvv(19:19)=='|')then
+ write(54,*)vvv(1:6), ' ', vvv(8:16), ' ', vvv(18:18)
+ go to 10
+ else if (vvv(20:20)=='|')then
+ write(54,*)vvv(1:6), ' ', vvv(8:16), ' ', vvv(18:19)
+ go to 10
+ else if (vvv(21:21)=='|')then
+ write(54,*)vvv(1:6), ' ', vvv(8:16), ' ', vvv(18:20)
+ go to 10
+ else if (vvv(22:22)=='|')then
+ write(54,*)vvv(1:6), ' ', vvv(8:16), ' ', vvv(18:21)
+ go to 10
+ else if (vvv(23:23)=='|')then
+ write(54,*)vvv(1:6), ' ', vvv(8:16), ' ', vvv(18:22)
+ go to 10
+ else if (vvv(25:25)=='|')then
+ write(54,*)vvv(1:6), ' ', vvv(8:16), ' ', vvv(18:24)
+ go to 10
+ endif
+ else if(vvv(18:18)=='|') then
+ if (vvv(20:20)=='|')then
+ write(54,*)vvv(1:6), ' ', vvv(8:17), ' ', vvv(19:19)
+ go to 10
+ else if (vvv(21:21)=='|')then
+ write(54,*)vvv(1:6), ' ', vvv(8:17), ' ', vvv(19:20)
+ go to 10
+ else if (vvv(22:22)=='|')then
+ write(54,*)vvv(1:6), ' ', vvv(8:17), ' ', vvv(19:21)
+ go to 10
+ else if (vvv(23:23)=='|')then
+ write(54,*)vvv(1:6), ' ', vvv(8:17), ' ', vvv(19:22)
+ go to 10
+ else if (vvv(24:24)=='|')then
+ write(54,*)vvv(1:6), ' ', vvv(8:17), ' ', vvv(19:23)
+ go to 10
+ else if (vvv(25:25)=='|')then
+ write(54,*)vvv(1:6), ' ', vvv(8:17), ' ', vvv(19:24)
+ go to 10
+ else if (vvv(26:26)=='|')then
+ write(54,*)vvv(1:6), ' ', vvv(8:17), ' ', vvv(19:25)
+ go to 10
+ endif
+ endif
+
+
+else if (vvv(9:9)=='|')then
+ if (vvv(17:17)=='|') then
+ if (vvv(24:24)=='|')then
+ write(54,*)vvv(1:8), ' ', vvv(10:16), ' ', vvv(18:23)
+ go to 10
+ else if (vvv(19:19)=='|')then
+ write(54,*)vvv(1:8), ' ', vvv(10:16), ' ', vvv(18:18)
+ go to 10
+ else if (vvv(20:20)=='|')then
+ write(54,*)vvv(1:8), ' ', vvv(10:16), ' ', vvv(18:19)
+ go to 10
+ else if (vvv(21:21)=='|')then
+ write(54,*)vvv(1:8), ' ', vvv(10:16), ' ', vvv(18:20)
+ go to 10
+ else if (vvv(22:22)=='|')then
+ write(54,*)vvv(1:8), ' ', vvv(10:16), ' ', vvv(18:21)
+ go to 10
+ else if (vvv(23:23)=='|')then
+ write(54,*)vvv(1:8), ' ', vvv(10:16), ' ', vvv(18:22)
+ go to 10
+ else if (vvv(25:25)=='|')then
+ write(54,*)vvv(1:8), ' ', vvv(10:16), ' ', vvv(18:24)
+ go to 10
+ endif
+ else if(vvv(19:19)=='|') then
+ if (vvv(21:21)=='|') then
+ write(54,*)vvv(1:8), ' ', vvv(10:18), ' ', vvv(20:20)
+ go to 10
+ else if (vvv(22:22)=='|') then
+ write(54,*)vvv(1:8), ' ', vvv(10:18), ' ', vvv(20:21)
+ go to 10
+ else if (vvv(23:23)=='|') then
+ write(54,*)vvv(1:8), ' ', vvv(10:18), ' ', vvv(20:22)
+ go to 10
+ else if (vvv(24:24)=='|') then
+ write(54,*)vvv(1:8), ' ', vvv(10:18), ' ', vvv(20:23)
+ go to 10
+ else if (vvv(25:25)=='|') then
+ write(54,*)vvv(1:8), ' ', vvv(10:18), ' ', vvv(20:24)
+ go to 10
+ else if (vvv(26:26)=='|') then
+ write(54,*)vvv(1:8), ' ', vvv(10:18), ' ', vvv(20:25)
+ go to 10
+ else if (vvv(27:27)=='|') then
+ write(54,*)vvv(1:8), ' ', vvv(10:18), ' ', vvv(20:26)
+ go to 10
+ endif
+ else if(vvv(20:20)=='|') then
+ if(vvv(22:22)=='|') then
+ write(54,*)vvv(1:8), ' ', vvv(10:19), ' ', vvv(21:21)
+ go to 10
+ else if (vvv(23:23)=='|')then
+ write(54,*)vvv(1:8), ' ', vvv(10:19), ' ', vvv(21:22)
+ go to 10
+ else if (vvv(24:24)=='|')then
+ write(54,*)vvv(1:8), ' ', vvv(10:19), ' ', vvv(21:23)
+ go to 10
+ else if (vvv(25:25)=='|')then
+ write(54,*)vvv(1:8), ' ', vvv(10:19), ' ', vvv(21:24)
+ go to 10
+ else if (vvv(26:26)=='|') then
+ write(54,*)vvv(1:8), ' ', vvv(10:19), ' ', vvv(21:25)
+ go to 10
+ else if (vvv(27:27)=='|') then
+ write(54,*)vvv(1:8), ' ', vvv(10:19), ' ', vvv(21:26)
+ go to 10
+ else if (vvv(28:28)=='|') then
+ write(54,*)vvv(1:8), ' ', vvv(10:19), ' ', vvv(21:27)
+ go to 10
+ endif
+ endif
+
+
+else if (vvv(10:10)=='|')then
+ if(vvv(18:18)=='|') then
+ if (vvv(20:20)=='|')then
+ write(54,*)vvv(1:9), ' ', vvv(11:17), ' ', vvv(19:19)
+ go to 10
+ else if (vvv(21:21)=='|')then
+ write(54,*)vvv(1:9), ' ', vvv(11:17), ' ', vvv(19:20)
+ go to 10
+ else if (vvv(22:22)=='|')then
+ write(54,*)vvv(1:9), ' ', vvv(11:17), ' ', vvv(19:21)
+ go to 10
+ else if (vvv(23:23)=='|')then
+ write(54,*)vvv(1:9), ' ', vvv(11:17), ' ', vvv(19:22)
+ go to 10
+ else if (vvv(24:24)=='|')then
+ write(54,*)vvv(1:9), ' ', vvv(11:17), ' ', vvv(19:23)
+ go to 10
+ else if (vvv(25:25)=='|')then
+ write(54,*)vvv(1:9), ' ', vvv(11:17), ' ', vvv(19:24)
+ go to 10
+ else if (vvv(26:26)=='|')then
+ write(54,*)vvv(1:9), ' ', vvv(11:17), ' ', vvv(19:25)
+ go to 10
+ endif
+ else if (vvv(20:20)=='|') then
+ if(vvv(22:22)=='|') then
+ write(54,*)vvv(1:9), ' ', vvv(11:19), ' ', vvv(21:21)
+ go to 10
+ else if (vvv(23:23)=='|')then
+ write(54,*)vvv(1:9), ' ', vvv(11:19), ' ', vvv(21:22)
+ go to 10
+ else if (vvv(24:24)=='|')then
+ write(54,*)vvv(1:9), ' ', vvv(11:19), ' ', vvv(21:23)
+ go to 10
+ else if (vvv(25:25)=='|')then
+ write(54,*)vvv(1:9), ' ', vvv(11:19), ' ', vvv(21:24)
+ go to 10
+ else if (vvv(26:26)=='|') then
+ write(54,*)vvv(1:9), ' ', vvv(11:19), ' ', vvv(21:25)
+ go to 10
+ else if (vvv(27:27)=='|') then
+ write(54,*)vvv(1:9), ' ', vvv(11:19), ' ', vvv(21:26)
+ go to 10
+ else if (vvv(28:28)=='|') then
+ write(54,*)vvv(1:9), ' ', vvv(11:19), ' ', vvv(21:27)
+ go to 10
+ endif
+ else if (vvv(21:21)=='|') then
+ if(vvv(23:23)=='|')then
+ write(54,*)vvv(1:9), ' ', vvv(11:20), ' ',vvv(22:22)
+ go to 10
+ else if(vvv(24:24)=='|')then
+ write(54,*)vvv(1:9), ' ', vvv(11:20), ' ',vvv(22:23)
+ go to 10
+ else if(vvv(25:25)=='|')then
+ write(54,*)vvv(1:9), ' ', vvv(11:20), ' ',vvv(22:24)
+ go to 10
+ else if(vvv(26:26)=='|')then
+ write(54,*)vvv(1:9), ' ', vvv(11:20), ' ',vvv(22:25)
+ go to 10
+ else if(vvv(27:27)=='|')then
+ write(54,*)vvv(1:9), ' ', vvv(11:20), ' ',vvv(22:26)
+ go to 10
+ else if(vvv(28:28)=='|')then
+ write(54,*)vvv(1:9), ' ', vvv(11:20), ' ',vvv(22:27)
+ go to 10
+ else if(vvv(29:29)=='|')then
+ write(54,*)vvv(1:9), ' ', vvv(11:20), ' ',vvv(22:28)
+ go to 10
+ endif
+ endif
+ go to 10
+endif
+
+11 close(53)
+close(54)
+
+end subroutine modifica_ascii
+
+
+
+
+
+!************************************************************************************************
+!Subroutine allocazione matrice sparsa
+!************************************************************************************************
+
+subroutine matrice_sparsa (nomefile,righe,colonne)
+!uso l'algoritmo CSR(Compressed )
+use dd
+implicit none
+character(40)::nomefile
+integer::i,j, cont, righe,colonne, nn, cont2
+real::a
+real,allocatable,dimension(:):: temp
+
+
+allocate(temp(colonne))
+allocate (iind(righe+1))
+open (unit=70, file=nomefile)
+cont=0
+
+do i=1,righe
+ read(70,*)(temp(j), j=1,colonne)
+ !write(6,*) temp
+ do j=1,colonne
+ if (temp(j)/=0.) then
+ cont=cont+1
+ endif
+ enddo
+enddo
+rewind(70)
+
+n_no_null=cont
+write(6,*)'There are', n_no_null, ' no-null cells in the Digital Terrain Model'
+
+allocate (val(n_no_null))
+allocate (jind(n_no_null))
+
+
+cont=0
+cont2=0
+do i=1,righe
+ cont2=cont+1
+ iind(i)=cont2
+ read (70,*) (temp(j),j=1,colonne)
+ do j=1,colonne
+ if (temp(j)/=0.) then
+ cont=cont+1
+ val(cont)=temp(j)
+ jind(cont)=j
+ endif
+ enddo
+enddo
+iind(righe+1)=n_no_null+1
+
+!write(6,*) "val=", val
+!write(6,*) "iind=", iind
+!write(6,*) "jind=", jind
+
+end subroutine matrice_sparsa
+
+
+
+
+
+!************************************************************************************************
+!Subroutine trova punto in matrice sparsa
+!************************************************************************************************
+
+subroutine trova_elemento(iriga,jcolonna)
+use dd
+integer::i,iriga,jcolonna, iniz, fine
+
+iniz=iind(iriga)
+fine=iind(iriga+1)-1
+
+do i=iniz,fine
+ if (jind(i)==jcolonna)then
+ output=val(i)
+ go to 144
+ else
+ output=0.
+ endif
+enddo
+
+144 i=0
+
+end subroutine trova_elemento
Added: grass-addons/raster/r.inund.fluv/find_main_channel.f90
===================================================================
--- grass-addons/raster/r.inund.fluv/find_main_channel.f90 (rev 0)
+++ grass-addons/raster/r.inund.fluv/find_main_channel.f90 2008-05-02 12:24:52 UTC (rev 31204)
@@ -0,0 +1,555 @@
+module dd
+integer,parameter::ferma_bisezione=15
+integer::n_no_null
+real, allocatable, dimension(:)::val, iind, jind
+real::output
+end module dd
+
+
+
+!*************************************************************************************************************************************
+!il programma in questione cerca dei limiti ragionevoli a destra e sinistra dell'alveo fluviale in questione
+!questi limiti sono quelli al di là dei quali ragionevolmente posso pensare che l'acqua è esondata
+!saranno poi i limiti al di là dei quali andrà a cercare il reticolo idrografico in base alla massima pendenza
+!COME AGISCE IL PROGRAMMA?
+!il programma noto il livello centennale ne calcola una porzione 30% (in maniera tale da escludere i difetti del dtm a fondo alveo)
+!quindi sale per step di mezzo metro e calcola i punti a dx o sx in cui l'acqua si scontra con il terreno (alveo)
+!si ferma al punto in cui a destra o a sinistra c'è una differenza sostanziale con il punto dello step prima
+!(paremetro da decidere = 1m, 2m)
+!tra i 2 passi cerco il punto del dtm con max pendenza "negativa" (quello è ¨il mio limite alveo)
+!INPUT
+!limiti della regione
+!matrice col DTM
+!profilo centennale
+!OUTPUT
+!vettoriale con i limiti dell'alveo
+!***************************************************************************************************************************************
+
+
+program trova_alveo
+use dd
+implicit none
+
+!****************************************************************************************************************************************
+!real, parameter:: res=2, par=2, step=0.5
+!il programma è pensato per una risoluzione del file raster 2x2 se la risoluzione fosse diversa MODIFICARE il parametro res
+!i valore 'par' e 'step' sono i valori su cui si basa il programma, modificare con ATTENZIONE
+!****************************************************************************************************************************************
+real::res, par, step
+real:: nord, sud, est, ovest
+integer::i, ii, iend, righe, colonne, nrighe, pixel_i, pixel_j, j, center_i, center_j, f
+integer::pixel_temp_i, pixel_temp_j, temp_x, temp_y
+!real, allocatable, dimension(:,:):: quota_terreno
+!character(2), allocatable, dimension(:,:):: limiti_alveo
+real(kind=8)::x, E1, N1, E1dopo, N1dopo, z, livello, E2, N2, dir_ortog, dist, E_temp, N_temp, coord_N, coord_E, quota, dist_temp, quota_max
+real(kind=8):: dist_destra, dist_sinistra, verifica_dist, pendenza, pendenza_max, quota_limite, quota_centro
+character(len=40)::xx, tab, tab1
+real,dimension(2)::quota_temp
+real, dimension(500):: E_ds, N_ds, E_sx, N_sx
+integer,allocatable,dimension(:)::i_ds, j_ds, i_sx, j_sx
+real,allocatable,dimension(:)::quota_ds, quota_sx, EE, NN
+character(40)::nomefile4,nomefile1, nomefile2, nomefile3
+open(unit=40, status='old', file='nomefile.txt')
+read(40,*) nomefile1, nomefile2, nomefile3, nomefile4
+close(40)
+
+! lettura parametri impostati dall'utente
+open(unit=41, status='old', file='parameter.txt')
+read(41,*)res, par, step
+close(41)
+
+!50=region2.txt
+!nomefile2=dtm2x2 (vd dopo)
+!51=quote_100
+!61=limiti_alveo_vect
+write(6,*)trim(nomefile1)
+open(unit=50, status='old', file=nomefile1)
+open(unit=51, status='old', file=nomefile4)
+open(unit=61, file=nomefile3)
+
+
+
+!leggo le dimensioni della regione utili per allocare le matrici in cui memorizzare carta di input e di output
+do i=1,13
+ read(50,*)xx
+ if (xx=='north:') then
+ backspace(50)
+ read(50,*)xx, nord
+ !write(60,'(a6,i)') trim(xx), nord
+ else if (xx=='south:') then
+ backspace(50)
+ read(50,*)xx, sud
+ !write(60,'(a6,i)')trim(xx), sud
+ else if (xx=='east:') then
+ backspace(50)
+ read(50,*)xx, est
+ !write(60,'(a5,i)') trim(xx), est
+ ! else if (xx=='nsres') then
+ ! backspace(50)
+ ! read(50,*)xx, res
+ else if (xx=='west:') then
+ backspace(50)
+ read(50,*)xx, ovest
+ !write(60,'(a5,i)')trim(xx), ovest
+ else if (xx=='rows:') then
+ backspace(50)
+ read(50,*)xx, righe
+ !write(60,'(a5,i)')trim(xx), righe
+ else if (xx=='cols:') then
+ backspace(50)
+ read(50,*)xx, colonne
+ !write(60,'(a5,i)')trim(xx), colonne
+ endif
+enddo
+write(6,*) nord, sud, est, ovest,res, righe, colonne
+rewind(50)
+close(50)
+
+!allocate(quota_terreno(righe,colonne))
+!allocate(limiti_alveo(righe,colonne))
+
+
+!leggo il DTM
+!read(52,*) ((quota_terreno(i,j), j=1,colonne), i=1,righe)
+!write(6,*) 'letto il DTM'
+!rewind(52)
+!close(52)
+
+call matrice_sparsa (nomefile2,righe,colonne)
+write(6,*) 'Read the Digital Terrain Model with Compressed Sparse Row (CSR) algorithm'
+
+
+!do j=1,colonne
+! do i=1,righe
+! limiti_alveo(i,j)='*'
+! enddo
+!enddo
+
+
+
+
+!guardo quant'è lungo il file che contiene il profilo
+do i=1,10000000
+ read(51,*,iostat=iend)xx
+ if(iend /= 0)exit
+enddo
+if(iend > 0)then
+ write(6,*)'Errore di lettura :',iend
+ stop
+endif
+nrighe=i-1
+write(6,*) 'Read',nrighe,' rows in file of water surface'
+rewind(51)
+
+
+allocate(i_ds(nrighe))
+allocate(j_ds(nrighe))
+allocate(quota_ds(nrighe))
+allocate(i_sx(nrighe))
+allocate(j_sx(nrighe))
+allocate(quota_sx(nrighe))
+allocate(EE(nrighe))
+allocate(NN(nrighe))
+
+
+
+do i=1,nrighe
+ i_ds(i)=0
+ j_ds(i)=0
+ quota_ds(i)=0
+ i_sx(i)=0
+ j_sx(i)=0
+ quota_sx(i)=0
+enddo
+
+i=0
+f=0
+
+!CICLO
+1000 i=i+1
+if (i==nrighe) then
+ go to 1001
+endif
+read(51,*)E1, N1
+!write(6,*)i, E1, N1
+read(51,*)E1dopo, N1dopo
+!write(6,*)i, E1dopo, N1dopo
+backspace(51)
+if (E1<ovest .or. E1>est .or. N1>nord .or. N1<sud) then
+ go to 1000
+endif
+f=f+1
+backspace(51)
+read(51,*) E1, N1, z
+!write(6,*)i, E1, N1, z
+!stop
+dir_ortog=-(E1dopo-E1)/(N1dopo-N1)
+if((N1dopo-N1)==0)then
+ dir_ortog=-99E10
+endif
+!write(6,*) dir_ortog
+!ora devo dire in quale pixel sono
+temp_x=int(E1)-ovest
+temp_y=nord-int(N1)
+center_i=1+int(temp_y/res)
+center_j=1+int(temp_x/res)
+call trova_elemento(center_i,center_j)
+quota_centro=output
+livello=z-quota_centro
+!cerco un secondo punto per costruire la retta
+E2=E1+10
+N2=N1+dir_ortog*(E2-E1)
+dist=sqrt((E2-E1)**2.+(N2-N1)**2.)
+!mi muovo sulla retta prima a ds e poi a sx
+
+
+!************************************************************************************************************************************************************
+!MI MUOVO A DESTRA
+j=0
+!elimino il 30%
+livello = 3./10.*livello
+2000 j=j+1
+livello=livello+step
+quota=livello+quota_centro
+!write(6,*)i, j, f, z, step, livello, quota, quota_terreno(center_i, center_j)
+if (quota>z)then
+ quota=z
+endif
+x=0
+2001 x=x+res
+coord_E=E1+x*(E2-E1)/dist
+coord_N=N1+x*(N2-N1)/dist
+temp_x=int(coord_E)-ovest
+temp_y=nord-int(coord_N)
+!ora devo dire in quale pixel sono
+pixel_i=1+int(temp_y/res)
+pixel_j=1+int(temp_x/res)
+!write(6,*) pixel_i, pixel_j, center_i, center_j
+if(pixel_i>righe .or. pixel_i<=0 .or. pixel_j<=0 .or. pixel_j>colonne)then
+ write(6,*) 'In' , f, "section of these region there isn't right boundary"
+ go to 2002
+endif
+call trova_elemento(pixel_i,pixel_j)
+if(output>=quota)then
+ E_ds(j)=1+int(coord_E/res)*res
+ N_ds(j)=1+int(coord_N/res)*res
+ !write(6,*)j, coord_E, E_ds(j), coord_N, N_ds(j)
+ if(j>1) then
+ dist_destra=sqrt((E_ds(j)-E_ds(j-1))**2.+(N_ds(j)-N_ds(j-1))**2)
+ !write(6,*)i,j,f,'sto andando a destra', dist_destra
+ !write(6,*) z, quota
+ if (dist_destra>par)then
+ !CERCO IL PUNTO IN CORRISPONDENZA DEL QUALE LA PENDENZA è MAX VERSO DESTRA '\ il punto ' è quello che considero come
+ !limite d'alveo oltre il quale calcolare il reticolo idrografico
+ !dovrò ancora verificare che tale pendenza sia effettivamente verso l'esterno alveo
+ !in caso contrario proseguo nel DTM fino a trovare un punto di pendenza verso l'esterno
+ E_temp=E_ds(j-1)
+ N_temp=N_ds(j-1)
+ pixel_temp_i=1+int((nord-N_temp)/res)
+ pixel_temp_j=int((E_temp-ovest)/res)+1
+ !write(6,*)pixel_temp_i, pixel_temp_j, righe, colonne, E_temp, N_temp
+ call trova_elemento(pixel_temp_i,pixel_temp_j)
+ quota_temp(1)=output
+ pendenza_max=100
+ verifica_dist=0
+20011 x=2
+ E_temp=E_temp+x*(E2-E1)/dist
+ N_temp=N_temp+x*(N2-N1)/dist
+ temp_x=int(E_temp)-ovest
+ temp_y=nord-int(N_temp)
+ pixel_temp_i=1+int(temp_y/res)
+ pixel_temp_j=1+int(temp_x/res)
+ if (pixel_temp_i<1 .or. pixel_temp_j<1 .or. pixel_temp_i>righe .or. pixel_temp_j>colonne)then
+ write(6,*) "In", f , "NO right boundary ~ the region is too small"
+ go to 2002
+ endif
+ call trova_elemento(pixel_temp_i,pixel_temp_j)
+ quota_temp(2)=output
+ if(quota_temp(2)==0)then
+ go to 20012
+ endif
+ pendenza=(quota_temp(2)-quota_temp(1))/x
+ if (pendenza<pendenza_max) then
+ i_ds(f)=pixel_temp_i
+ j_ds(f)=pixel_temp_j
+ EE(f)=E_temp
+ NN(f)=N_temp
+ quota_ds(f)=quota_temp(2)
+ pendenza_max=pendenza
+ endif
+ quota_temp(1)=quota_temp(2)
+ verifica_dist=verifica_dist+x
+ if(verifica_dist<dist_destra)then
+ go to 20011
+ else
+ if(pendenza_max<=0) then
+ write(61,*) EE(f), NN(f), quota_ds(f), z, f
+ go to 2002
+ else
+ go to 20011
+ endif
+ endif
+ else
+ if(quota==z)then
+ write(6,*) 'In', f, "section NO right boundary ~ SEE THE PARAMETER"
+ go to 2002
+ endif
+ go to 2000
+ endif
+ else
+ go to 2000
+ endif
+else
+ go to 2001
+endif
+
+
+20012 if(pendenza_max<=0) then
+ write(61,*) EE(f), NN(f), quota_ds(f), z, f
+ go to 2002
+else
+ write(6,*) "In", f, "section the boundary of river is out of DTM"
+ go to 2002
+endif
+
+
+
+
+!************************************************************************************************************************************************************
+!MI MUOVO A SINISTRA
+
+2002 livello=z-quota_centro
+j=0
+!elimino il 30%
+livello = 3/10*livello
+2003 j=j+1
+livello=livello+step
+quota=livello+quota_centro
+!write(6,*)i, j, f, z, step, livello, quota, quota_terreno(center_i, center_j)
+!stop
+if (quota>z)then
+ quota=z
+endif
+x=-2
+2004 x=x-2
+coord_E=E1+x*(E2-E1)/dist
+coord_N=N1+x*(N2-N1)/dist
+temp_x=int(coord_E)-ovest
+temp_y=nord-int(coord_N)
+!ora devo dire in quale pixel sono
+pixel_i=1+int(temp_y/res)
+pixel_j=1+int(temp_x/res)
+!write(6,*) pixel_i, pixel_j, center_i, center_j
+if(pixel_i>righe .or. pixel_i<=0 .or. pixel_j<=0 .or. pixel_j>colonne)then
+ write(6,*)'In' , f, "section of these region there isn't left boundary"
+ go to 1000
+endif
+call trova_elemento(pixel_i,pixel_j)
+if(output>=quota)then
+ E_sx(j)=1+int(coord_E/res)*res
+ N_sx(j)=1+int(coord_N/res)*res
+ !write(6,*) coord_E, E_ds(j), coord_N, N_ds(j)
+ if(j>1) then
+ dist_sinistra=sqrt((E_sx(j)-E_sx(j-1))**2.+(N_sx(j)-N_sx(j-1))**2)
+ !write(6,*) i,j,f, 'sto andando a sinistra',dist_sinistra
+ !write(6,*) z, quota
+ if (dist_sinistra>par)then
+ !CERCO IL PUNTO IN CORRISPONDENZA DEL QUALE LA PENDENZA è MAX VERSO SINISTRA /' il punto ' è quello che considero come
+ !limite d'alveo oltre il quale calcolare il reticolo idrografico
+ !dovrò ancora verificare che tale pendenza sia effettivamente verso l'esterno alveo
+ !in caso contrario proseguo nel DTM fino a trovare un punto di pendenza verso l'esterno
+ E_temp=E_sx(j-1)
+ N_temp=N_sx(j-1)
+ pixel_temp_i=1+int((nord-N_temp)/res)
+ pixel_temp_j=int((E_temp-ovest)/res)+1
+ call trova_elemento(pixel_temp_i,pixel_temp_j)
+ quota_temp(1)=output
+ pendenza_max=100
+ verifica_dist=0
+20041 x=-res
+ E_temp=E_temp+x*(E2-E1)/dist
+ N_temp=N_temp+x*(N2-N1)/dist
+ temp_x=int(E_temp)-ovest
+ temp_y=nord-int(N_temp)
+ pixel_temp_i=1+int(temp_y/res)
+ pixel_temp_j=1+int(temp_x/res)
+ if (pixel_temp_i<1 .or. pixel_temp_j<1 .or. pixel_temp_i>righe .or. pixel_temp_j>colonne)then
+ write(6,*) "In", f , "NO left boundary ~ the region is too small"
+ go to 1000
+ endif
+ call trova_elemento(pixel_temp_i,pixel_temp_j)
+ quota_temp(2)=output
+ if(quota_temp(2)==0)then
+ go to 20042
+ endif
+ pendenza=(quota_temp(2)-quota_temp(1))/x
+ if (pendenza<pendenza_max) then
+ i_sx(f)=pixel_temp_i
+ j_sx(f)=pixel_temp_j
+ EE(f)=E_temp
+ NN(f)=N_temp
+ quota_sx(f)=quota_temp(2)
+ pendenza_max=pendenza
+ endif
+ quota_temp(1)=quota_temp(2)
+ verifica_dist=verifica_dist+x
+ if(abs(verifica_dist)<dist_destra)then
+ go to 20041
+ else
+ if(pendenza_max<=0) then
+ write(61,*) EE(f), NN(f), quota_sx(f), z, f
+ go to 1000
+ else
+ go to 20041
+ endif
+ endif
+ else
+ if(quota==z)then
+ write(6,*)'In', f, "section NO left boundary ~ SEE THE PARAMETER"
+ go to 1000
+ endif
+ go to 2003
+ endif
+ else
+ go to 2003
+ endif
+else
+ go to 2004
+endif
+
+20042 if(pendenza_max<=0) then
+ write(61,*) EE(f), NN(f), quota_sx(f), z, f
+ go to 1000
+else
+ write(6,*)"In", f, "section the boundary of river is out of DTM"
+ go to 1000
+endif
+
+
+
+
+1001 write(6,*) "End of water surface file"
+
+
+
+
+write(6,*) 'Writing the boundary profile'
+
+!write(60,*) ((limiti_alveo(i,j), j=1,colonne), i=1,righe)
+
+deallocate(i_ds)
+deallocate(j_ds)
+deallocate(quota_ds)
+deallocate(i_sx)
+deallocate(j_sx)
+deallocate(quota_sx)
+deallocate(EE)
+deallocate(NN)
+
+write(6,*) 'END of Fortran Code ~ OK'
+end program trova_alveo
+
+
+
+
+
+
+
+
+
+!************************************************************************************************
+!Subroutine allocazione matrice sparsa
+!************************************************************************************************
+
+subroutine matrice_sparsa (nomefile,righe,colonne)
+!uso l'algoritmo CSR(Compressed )
+use dd
+implicit none
+character(40)::nomefile
+integer::i,j, cont, righe,colonne, nn, cont2
+real::a
+real,allocatable,dimension(:):: temp
+
+
+allocate(temp(colonne))
+allocate (iind(righe+1))
+
+open (unit=70, file=nomefile)
+cont=0
+do i=1,righe
+ read(70,*)(temp(j), j=1,colonne)
+ !write(6,*) temp
+ do j=1,colonne
+ if (temp(j)/=0.) then
+ cont=cont+1
+ endif
+ enddo
+enddo
+rewind(70)
+
+n_no_null=cont
+write(6,*)'There are', n_no_null, 'no-null cells in Digital Terrain Model '
+
+allocate (val(n_no_null))
+allocate (jind(n_no_null))
+
+
+cont=0
+cont2=0
+do i=1,righe
+ cont2=cont+1
+ iind(i)=cont2
+ read (70,*) (temp(j),j=1,colonne)
+ do j=1,colonne
+ if (temp(j)/=0.)then
+ cont=cont+1
+ val(cont)=temp(j)
+ jind(cont)=j
+ endif
+ enddo
+enddo
+iind(righe+1)=n_no_null+1
+
+rewind(70)
+close(70)
+
+!write(6,*) "val=", val
+!write(6,*) "iind=", iind
+!write(6,*) "jind=", jind
+!stop
+
+end subroutine matrice_sparsa
+
+
+
+
+
+!************************************************************************************************
+!Subroutine trova punto in matrice sparsa
+!************************************************************************************************
+
+subroutine trova_elemento(iriga,jcolonna)
+use dd
+implicit none
+integer::i,iriga,jcolonna, iniz, fine
+
+if (iind(iriga)==iind(iriga+1))then
+ output=0.
+ go to 144
+endif
+iniz=iind(iriga)
+fine=iind(iriga+1)-1
+
+do i=iniz,fine
+ if (jind(i)==jcolonna)then
+ output=val(i)
+ go to 144
+ else
+ output=0.
+ endif
+enddo
+
+144 i=0
+!write(6,*) output,iriga, jcolonna, iniz, fine
+
+end subroutine trova_elemento
Added: grass-addons/raster/r.inund.fluv/install.sh
===================================================================
--- grass-addons/raster/r.inund.fluv/install.sh (rev 0)
+++ grass-addons/raster/r.inund.fluv/install.sh 2008-05-02 12:24:52 UTC (rev 31204)
@@ -0,0 +1,62 @@
+#!/bin/sh
+
+echo "You must be root or have the permissions for installing a grass-script"
+echo " "
+echo "please, write your GRASS directory"
+echo " e.g. /usr/grass-xxxx (look also in /usr/local or /usr/lib)"
+read grass_directory
+
+if [ -e $grass_directory/fortran_code ]; then
+ directory=$grass_directory'/fortran_code'
+else
+ echo "The fortran code will be put in $grass_directory/fortran_code"
+ directory=$grass_directory'/fortran_code'
+ mkdir $directory
+fi
+
+initial_directory=`pwd`
+
+cp find_main_channel.f90 $directory
+cp clean_inundation.f90 $directory
+cp 2d_path.f90 $directory
+cp correction_from_path.f90 $directory
+
+cd $directory
+
+
+# if you have gnu fortran compiler (gcc-gfortran.i386), use the following rows:
+
+gfortran -O1 -o find_main_channel.exe find_main_channel.f90
+gfortran -O1 -o clean_inundation.exe clean_inundation.f90
+gfortran -O1 -o 2d_path.exe 2d_path.f90
+gfortran -O1 -o correction_from_path.exe correction_from_path.f90
+
+# if you have a intel fortran compiler, use the following rows
+#ifort -O3 -xW -o find_main_channel.exe find_main_channel.f90
+#ifort -O3 -xW -o clean_inundation.exe clean_inundation.f90
+#ifort -O3 -xW -o 2d_path.exe 2d_path.f90
+#ifort -O3 -xW -o correction_from_path.exe correction_from_path.f90
+
+rm find_main_channel.f90
+rm clean_inundation.f90
+rm 2d_path.f90
+rm correction_from_path.f90
+
+cd $initial_directory
+
+cp r.inund.fluv $grass_directory/scripts
+chmod 777 $grass_directory/scripts/r.inund.fluv
+
+if [ -e $grass_directory/docs ]; then
+ cp r.inund.fluv.html $grass_directory/docs/html/
+ else
+ mkdir $grass_directory/docs
+ mkdir $grass_directory/docs/html
+ cp r.inund.fluv.html $grass_directory/docs/html/
+fi
+
+echo "OK now you've successfull installed "
+
+exit
+
+
Added: grass-addons/raster/r.inund.fluv/r.inund.fluv
===================================================================
--- grass-addons/raster/r.inund.fluv/r.inund.fluv (rev 0)
+++ grass-addons/raster/r.inund.fluv/r.inund.fluv 2008-05-02 12:24:52 UTC (rev 31204)
@@ -0,0 +1,507 @@
+#!/bin/sh
+############################################################################
+#
+# MODULE: r.inund.fluv v.1 for GRASS 6.3 (april 2008)
+# AUTHOR(S): Roberto Marzocchi, Bianca Federici, Domenico Sguerso
+# PURPOSE: Creates a fluvial inundation map given an high-resolution dtm and a water surface profile
+# COPYRIGHT: (C) 2008 by the GRASS Development Team
+#
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+#############################################################################
+
+#%Module
+#% description: Creates a fluvial inundation map given an high-resolution dtm and a water surface profile
+#% keywords: Automatic procedure to compute a fluvial inundation map
+#%End
+#%option
+#% key: DTM
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Input DTM raster map
+#% required : yes
+#%end
+#%option
+#% key: W_S_PROFILE
+#% type: string
+#% gisprompt: old_file,file,input
+#% description: Input ASCII file of the water surface profile
+#% required : yes
+#%end
+#%option
+#% key: RIVER
+#% type: string
+#% gisprompt: old,vector,vector
+#% description: Input vector line map of river-axis
+#% required : yes
+#%end
+#%option
+#% key: FLOODING_MAP
+#% type: string
+#% gisprompt: new,cell,raster
+#% description: Output: name of flooding map
+#% required : yes
+#%end
+#%option
+#% key: DOUBT_MAP
+#% type: string
+#% gisprompt: new,cell,raster
+#% description: Output: name of doubful surface areas
+#% required : no
+#%end
+#%option
+#%guisection:Optional parameters & maps
+#% key: PROFILE_T100
+#% type: string
+#% gisprompt: old_file,file,input
+#% description: Input ASCII file with water-depht for return period T > 100 years
+#% required : no
+#%end
+#%option
+#%guisection:Optional parameters & maps
+ #% key: delta_y
+#% type: double
+#% gisprompt: new
+#% description: Input delta_y to find the boundaries of the main channel [default value 0.5 m]
+#% required : no
+#%end
+#%option
+#%guisection:Optional parameters & maps
+#% key: delta_x
+#% type: double
+#% gisprompt: new
+#% description: Input delta_x to find the boundaries of the main channel [default value 3.5 m]
+#% required : no
+#%end
+#%option
+#%guisection:Optional parameters & maps
+#% key: river_boundary
+#% type: string
+#% gisprompt: new,vector,vector
+#% description: Output vector boundaries of the main channel
+#% required : no
+#%end
+#%option
+#%guisection:Optional parameters & maps
+#% key: boundary_type
+#% type: string
+#% description: Format of vector boundaries of the main channel
+#% options: line,points
+#% answer: points
+#% required: no
+#%end
+#%option
+#%guisection:Optional parameters & maps
+#% key: res_B
+#% type: integer
+#% gisprompt: new,cell,raster
+#% description: Input value: resolution B [default value 10 m]
+#% required : no
+#%end
+#%option
+#%guisection:Optional parameters & maps
+#% key: res_C
+#% type: integer
+#% gisprompt: new,cell,raster
+#% description: Input value: resolution C [default value 20 m]
+#% required : no
+#%end
+#%option
+#%guisection:Optional parameters & maps
+#% key: MAP_W_S_PROFILE
+#% type: string
+#% gisprompt: new,vector,vector
+#% description: Output vector point map of the water surface profile
+#% required : no
+#%end
+
+if [ -z "$GISBASE" ]
+then
+ echo ""
+ echo "You must be in GRASS GIS to run this program"
+ echo ""
+ exit 1
+fi
+
+if [ "$1" != "@ARGS_PARSED@" ] ; then
+ exec g.parser "$0" "$@"
+fi
+
+$GRASS_VERBOSE
+PROG=`basename $0`
+#request control & test (only for imput map)
+if [ -z "$GIS_OPT_DTM" ] ; then
+ g.message -e "ERROR: DTM not specified"
+ exit 1
+fi
+#g.findfile element=cell file="$GIS_OPT_DTM" > /dev/null # test if a map exists
+#if [ $? -eq 0 ]; then
+# echo "Imput map not found" 1>&2
+# exit 1
+#fi
+
+g.region save=dtm_res2 --overwrite
+g.region rast="$GIS_OPT_DTM" save=verifica
+g.region -p > .resolution.txt
+g.region region=dtm_res2
+exec 6<&0
+exec < .resolution.txt
+read a1 b1
+read a2 b2
+read a3 b3
+read a4 b4
+read a5 b5
+read a6 b6
+read a7 b7
+read a8 b8
+read a9 b9
+exec 0<&6 6<&-
+res=$b9
+if [ "$b9" -gt 5 ]
+then
+ g.message -w "WARNING - the DTM resolution it'snt good for this application"
+else
+ g.message "The DTM resolution is $res x $res meters"
+fi
+
+
+rm .resolution.txt
+
+
+if [ -z "$GIS_OPT_W_S_PROFILE" ] ; then
+ g.message -e "ERROR: file with depht not specified"
+ exit 1
+fi
+
+
+
+if [ -z "$GIS_OPT_RIVER" ] ; then
+ g.message -e "ERROR: vector with river-axis not specified"
+ exit 1
+fi
+
+
+if [ -z "$GIS_OPT_FLOODING_MAP" ] ; then
+ g.message -e "ERROR: name of flooding map not specified"
+ exit 1
+fi
+
+
+#optional control: assignement default value
+if [ -n "$GIS_OPT_res_B" ]
+then
+ res10="$GIS_OPT_res_B"
+else
+ res10=10
+fi
+
+if [ -n "$GIS_OPT_res_C" ]
+then
+ res20="$GIS_OPT_res_C"
+else
+ res20=20
+fi
+
+if [ -n "$GIS_PROFILE_T100" ]
+then
+ quote100="$GIS_OPT_PROFILE_T100"
+else
+ quote100="$GIS_OPT_W_S_PROFILE"
+fi
+
+if [ -n "$GIS_OPT_delta_x" ]
+then
+ deltax="$GIS_OPT_delta_x"
+else
+ deltax=$(echo "$res * 1.5" | bc)
+fi
+
+if [ -n "$GIS_OPT_delta_y" ]
+then
+ deltay="$GIS_OPT_delta_y"
+else
+ deltay="0.5"
+fi
+
+g.remove rast=MASK --quiet
+cartella_nascosta=".temp_grass_script"
+# do not use any space in name of cartella_nascosta
+
+
+rm -r -f ~/$cartella_nascosta
+mkdir ~/$cartella_nascosta
+cp "$GISBASE/fortran_code"/find_main_channel.exe ~/$cartella_nascosta
+cp "$GISBASE/fortran_code"/clean_inundation.exe ~/$cartella_nascosta
+cp "$GISBASE/fortran_code"/2d_path.exe ~/$cartella_nascosta
+cp "$GISBASE/fortran_code"/correction_from_path.exe ~/$cartella_nascosta
+
+#variables
+dtm="$GIS_OPT_DTM"
+cp "$GIS_OPT_W_S_PROFILE" ~/$cartella_nascosta/profilo
+v.in.ascii -z input="$GIS_OPT_W_S_PROFILE" output=profilo format=point fs=" " skip=0 x=1 y=2 z=3 cat=0 --overwrite --quiet
+centerline="$GIS_OPT_RIVER"
+inondazione="$GIS_OPT_FLOODING_MAP"
+
+
+#calculation of the "boundaries" of the main channel
+g.region nsres=$res ewres=$res save=dtm_res2 --overwrite
+g.region -p >> ~/$cartella_nascosta/"region2.txt"
+r.out.ascii -h input=$dtm output=$cartella_nascosta/"dtm2x2" null=0 --quiet
+
+#fortran code run
+cd
+cp "$quote100" ~/$cartella_nascosta/quote_100
+cd $cartella_nascosta
+echo $res $deltax $deltay > parameter.txt
+echo "'"region2.txt"'" "'"dtm2x2"'" "'"limiti_alveo_vect"'" "'"quote_100"'" > nomefile.txt
+g.message "Fortran code that automatically find bed river"
+./find_main_channel.exe
+rm nomefile.txt
+rm parameter.txt
+#TEST FORTRAN CODE RUN
+if [ -z "limiti_alveo_vect" ]
+then
+ g.message -e "ERROR IN FORTRAN CODE - see fortran error message and try understand the problem"
+ exit 1
+fi
+cd
+v.in.ascii --o input=$cartella_nascosta/limiti_alveo_vect output=limiti_alveo format=point fs=' ' skip=0 x=1 y=2 cat=0 --quiet
+
+
+#export of 20x20 dtm, region20.txt and region10.txt
+g.region nsres=$res20 ewres=$res20 save=dtm_res20 --overwrite
+g.region -p >> ~/$cartella_nascosta/region20.txt
+r.out.ascii -h input=$dtm output=$cartella_nascosta/dtm20x20 null=0 --quiet
+g.region nsres=$res10 ewres=$res10 save=dtm_res10 --overwrite
+g.region -p >> ~/$cartella_nascosta/region10.txt
+
+
+g.message "STEP 1"
+###################################################################################################################
+# STEP 1 - Thiessen interpolation
+###################################################################################################################
+# don't use low resolution [default 10m]
+r.mask input=$dtm maskcats=* --quiet
+g.region region=dtm_res10
+v.surf.idw input=profilo output=quote_punti_adiacenti npoints=1 layer=1 column=dbl_3 --quiet
+g.region region=dtm_res2
+r.mapcalculator amap=quote_punti_adiacenti bmap=$dtm formula="if(A>=B)" outfile=inondazione_step1 help=- --quiet
+
+
+g.message "STEP 2"
+###################################################################################################################
+# STEP 2 - remove "lakes"
+###################################################################################################################
+r.to.vect -s input=inondazione_step1 output=inondazione_step1 feature=area --overwrite --quiet
+v.select ainput=inondazione_step1 atype=line,area alayer=1 binput=$centerline btype=line blayer=1 output=inondazione_step2 operator=overlap --overwrite --quiet
+v.to.rast input=inondazione_step2 output=inondazione_step2 use=val layer=1 value=1 rows=4096 --overwrite --quiet
+
+
+g.message "STEP 3"
+###################################################################################################################
+# STEP 3 - looking at the minimum path from wet pixel to the channel axis
+# if the dtm elevation is higher than the water elevation in the pixel ==> pixel becomes dry
+# if the dtm elevation is always lower than the water elevation in the pixel ==> pixel remains wet
+###################################################################################################################
+# export data
+#20x20
+g.region region=dtm_res20
+r.out.ascii -h input=inondazione_step2 output=$cartella_nascosta/inondazione_step2 null=0 --quiet
+r.out.ascii -h input=quote_punti_adiacenti output=$cartella_nascosta/quote_punti_adiacenti dp=2 null=0 --quiet
+
+# FORTRAN code
+g.message "CODE FORTRAN OF STEP 3"
+cd
+cd $cartella_nascosta
+echo "'"inondazione_step2"'" "'"profilo"'" "'"inondazione_nuova.txt"'" "'"quote_punti_adiacenti"'" "'"region20.txt"'" "'"region2.txt"'" "'"dtm2x2"'" > nomefile.txt
+./clean_inundation.exe
+rm nomefile.txt
+
+#TEST FORTRAN CODE RUN
+if [ -z "inondazione_nuova.txt" ]
+then
+ g.message -e "ERROR IN FORTRAN CODE - see fortran error message and try understand the problem"
+ exit 1
+fi
+cd
+
+# output file
+r.in.ascii input=$cartella_nascosta/inondazione_nuova.txt output=inondazione_step3 'mult=1.0 or read from header' nv=0 --overwrite --quiet
+
+
+g.message "STEP 4"
+###################################################################################################################
+# STEP 4 - calcolation of the 2d_path outside the main channel
+###################################################################################################################
+r.out.ascii -h input=inondazione_step3 output=$cartella_nascosta/inondazione_step3 null=0 --quiet
+
+g.message "2D MODEL OF DIFFUSION OF FLOODING INONDATION OUTSIDE THE MAIN CHANEL"
+cd
+cd $cartella_nascosta
+echo "'"profilo"'" "'"inondazione_step3"'" "'"reticolo"'" "'"region20.txt"'" "'"dtm20x20"'" "'"limiti_alveo_vect"'" "'"region2.txt"'" "'"dtm2x2"'" "'"bordi_alveo"'" > nomefile.txt
+./2d_path.exe
+rm nomefile.txt
+
+#TEST FORTRAN CODE RUN
+if [ -z "reticolo" ]
+then
+ g.message -e "ERROR IN FORTRAN CODE - see fortran error message and try understand the problem"
+ exit 1
+fi
+cd
+
+
+#r.in.ascii input=$cartella_nascosta/bordi_alveo output=bordi_alveo 'mult=1.0 or read from header' nv=0 --overwrite
+#r.thin input=bordi_alveo output=bordi_alveo_thin iterations=200 --overwrite
+#r.to.vect input=bordi_alveo_thin output=bordi_alveo feature=line --overwrite
+v.in.ascii -n input=$cartella_nascosta/bordi_alveo output=bordi_alveo format=standard fs=' ' skip=0 x=1 y=2 z=0 cat=0 --quiet
+r.in.ascii input=$cartella_nascosta/reticolo output=reticolo 'mult=1.0 or read from header' nv=* --overwrite --quiet
+r.thin input=reticolo output=reticolo_thin iterations=200 --overwrite --quiet
+r.to.vect input=reticolo_thin output=reticolo feature=line --overwrite --quiet
+g.remove rast=reticolo_thin --quiet
+
+
+# 10x10 resolution
+---> 1 pò più lento ma più preciso!
+g.region region=dtm_res10
+
+# difference between maps of step1 and step2
+r.mapcalculator amap=inondazione_step3 formula='isnull(A)' outfile=inondazione_step3_isnull help=- --quiet
+r.mapcalculator amap=inondazione_step3_isnull formula='if(A,0,1)' outfile=inondazione_step3 help=- --quiet
+r.mapcalculator amap=inondazione_step2 bmap=inondazione_step3 formula='A-B' outfile=diff_inond help=- --quiet
+r.mapcalculator amap=diff_inond formula='if(A,1,0,0)' outfile=diff_inond help=- --quiet #remove the -1 value resulted from A-B subtraction
+r.null map=diff_inond setnull=0 --quiet
+r.to.vect input=diff_inond output=diff_inond1 feature=area --overwrite --quiet
+#to use the v.select command the line must have a layer associated ==> export of vector map in dxf format and then import the dxf file
+v.out.dxf input=reticolo output=$cartella_nascosta/reticolo.dxf --quiet
+v.in.dxf -1 input=$cartella_nascosta/reticolo.dxf output=reticolo --overwrite --quiet
+
+#only 2d_path overlapped the surfaces to look
+v.select ainput=reticolo atype=line alayer=1 binput=diff_inond1 btype=area blayer=1 output=reticolo_interno operator=overlap --overwrite --quiet
+
+# creation of points by 2d path outside the main channel
+v.to.points input=reticolo_interno type=point,line,boundary,centroid output=punti_reticolo_interno llayer=1 dmax=30 --overwrite --quiet
+
+# convert 2d points to 3d points with the elevation of the water surface profile at the starting point of the path in the channel axis (or with the elevation of the nearest water surface profile)
+g.region region=dtm_res20 --quiet
+v.drape input=punti_reticolo_interno type=point,centroid,line,boundary,face,kernel rast=reticolo method=nearest output=punti3d_reticolo --overwrite --quiet
+#v.drape input=punti_reticolo_interno type=point,centroid,line,boundary,face,kernel rast=quote_punti_adiacenti method=nearest output=punti3d_reticolo --overwrite --quiet
+
+
+g.region region=dtm_res10
+
+#export of 3d points
+v.out.ascii input=punti3d_reticolo output=$cartella_nascosta/punti3d_reticolo format=point dp=2 --quiet
+r.out.ascii -h input=diff_inond output=$cartella_nascosta/diff_inond null=0 --quiet
+
+
+
+g.message "FORTRAN CODE OF STEP 4"
+cd
+cd $cartella_nascosta
+echo "'"diff_inond"'" "'"punti3d_reticolo"'" "'"correzione"'" "'"region10.txt"'" "'"dtm2x2"'" "'"region2.txt"'" > nomefile.txt
+./correction_from_path.exe
+rm nomefile.txt
+
+#TEST FORTRAN CODE RUN
+if [ -z "correzione" ]
+then
+ g.message -e "ERROR IN FORTRAN CODE - see fortran error message and try understand the problem"
+ exit 1
+fi
+cd
+
+
+r.in.ascii input=$cartella_nascosta/correzione output=inondazione_step5 'mult=1.0 or read from header' nv=0 --overwrite --quiet
+
+# replace null cells with 0 cells
+g.region region=dtm_res2
+r.mapcalculator amap=inondazione_step5 formula='isnull(A)' outfile=inondazione_step5_isnull help=- --quiet
+r.mapcalculator amap=inondazione_step5_isnull formula='if(A,0,1)' outfile=inondazione_step5 help=- --quiet
+r.mapcalculator amap=inondazione_step5 bmap=inondazione_step3 formula='A+B' outfile=inondazione help=- --quiet
+
+
+g.message "STEP 5"
+###################################################################################################################
+# STEP 5 - EXPORT OUTPUT FILES
+###################################################################################################################
+g.rename rast=inondazione,inondazione1 --overwrite --quiet
+r.null map=inondazione1 setnull=0 --quiet
+r.to.vect -s input=inondazione1 output=inondazione1 feature=area --overwrite --quiet
+v.patch input=centerline,punti3d_reticolo output=overlap --overwrite --quiet
+v.select ainput=inondazione1 atype=area alayer=1 binput=overlap btype=line blayer=1 output=inondazione operator=overlap --overwrite --quiet
+v.to.rast input=inondazione output=$inondazione use=val layer=1 value=1 rows=4096 --overwrite --quiet
+
+#COLOURS OF INUNDATION RASTER MAP
+rm -f ".colore_inondazione.file"
+echo '1 blue' >> ".colore_inondazione.file"
+echo '0 grey' >> ".colore_inondazione.file"
+echo 'nv white' >> ".colore_inondazione.file"
+echo 'end' >> ".colore_inondazione.file"
+cat ".colore_inondazione.file" | r.colors map=$inondazione color=rules
+rm ".colore_inondazione.file"
+
+
+
+if [ -n "$GIS_OPT_DOUBT_MAP" ]
+then
+ # It's possible to create a map with doubful surface areas
+ dubbio="$GIS_OPT_DOUBT_MAP"
+ r.mapcalculator amap=$inondazione formula='isnull(A)' outfile=temp help=- --quiet
+ r.mapcalculator amap=temp formula='if(A,0,1)' outfile=temp help=- --quiet
+ r.mapcalculator amap=inondazione_step2 bmap=temp formula='A-B' outfile=$dubbio help=- --quiet
+ r.null map=$dubbio setnull=0 --quiet
+ #COLOURS OF DOUBFUL RASTER MAP
+ rm -f .colore_dubbio.file
+ echo '1 red' >> .colore_dubbio.file
+ echo '0 grey' >> .colore_dubbio.file
+ echo 'nv white' >> .colore_dubbio.file
+ echo 'end' >> .colore_dubbio.file
+ cat .colore_dubbio.file | r.colors map=$dubbio color=rules
+ rm .colore_dubbio.file
+ g.remove rast=temp
+fi
+
+
+g.remove rast=MASK --quiet
+g.remove rast=inondazione_step1,inondazione_step2,inondazione_step3,reticolo,quote_punti_adiacenti --quiet
+g.remove rast=inondazione_step3_isnull,diff_inond,inondazione_step5,inondazione_step5_isnull,inondazione1 --quiet
+g.remove vect=inondazione_step1,inondazione_step2,reticolo --quiet
+g.remove vect=diff_inond1,reticolo_interno,inondazione1,overlap --quiet
+
+if [ "$inondazione" == "inondazione" ]
+then
+ g.message " "
+else
+ g.remove vect=inondazione --quiet
+fi
+
+if [ -n "$GIS_OPT_MAP_W_S_PROFILE" ]
+then
+ g.rename vect=profilo,"$GIS_OPT_MAP_W_S_PROFILE" --quiet
+else
+ g.remove vect=profilo --quiet
+fi
+
+
+if [ -n "$GIS_OPT_river_boundary" ]
+then
+ if [ "$GIS_OPT_boundary_type" = "line" ]
+ then
+ g.rename vect=bordi_alveo,"$GIS_OPT_river_boundary"
+ g.remove vect=limiti_alveo --quiet
+ else
+ g.rename vect=limiti_alveo,"$GIS_OPT_river_boundary"
+ g.remove vect=bordi_alveo --quiet
+ fi
+else
+ g.remove vect=bordi_alveo,limiti_alveo --quiet
+fi
+
+rm -r ~/$cartella_nascosta
+
+
+
+exit 0
Added: grass-addons/raster/r.inund.fluv/r.inund.fluv.html
===================================================================
--- grass-addons/raster/r.inund.fluv/r.inund.fluv.html (rev 0)
+++ grass-addons/raster/r.inund.fluv/r.inund.fluv.html 2008-05-02 12:24:52 UTC (rev 31204)
@@ -0,0 +1,148 @@
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html>
+<head>
+<title>GRASS GIS: r.inund.fluv</title>
+<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
+<link rel="stylesheet" href="grassdocs.css" type="text/css">
+</head>
+<body bgcolor="white">
+
+<img src="grass_logo.png" alt="GRASS logo"><hr align=center size=6 noshade>
+
+<h2>NAME</h2>
+<em><b>r.inund.fluv</b></em> - Creates a fluvial inundation map given an high-resolution dtm and a water surface profile
+<h2>KEYWORDS</h2>
+Automatic procedure to compute a fluvial inundation map
+<h2>SYNOPSIS</h2>
+<b>r.inund.fluv</b><br>
+<b>r.inund.fluv help</b><br>
+<b>r.inund.fluv</b> <b>DTM</b>=<em>string</em> <b>W_S_PROFILE</b>=<em>string</em> <b>RIVER</b>=<em>string</em> <b>FLOODING_MAP</b>=<em>string</em> [<b>DOUBT_MAP</b>=<em>string</em>] [<b>PROFILE_T100</b>=<em>string</em>] [<b>(null)</b>=<em>float</em>] [<b>delta_x</b>=<em>float</em>] [<b>river_boundary</b>=<em>string</em>] [<b>boundary_type</b>=<em>string</em>] [<b>res_B</b>=<em>integer</em>] [<b>res_C</b>=<em>integer</em>] [<b>MAP_W_S_PROFILE</b>=<em>string</em>] [--<b>overwrite</b>] [--<b>verbose</b>] [--<b>quiet</b>]
+
+<h3>Flags:</h3>
+<DL>
+<DT><b>--overwrite</b></DT>
+<DD>Allow output files to overwrite existing files</DD>
+<DT><b>--verbose</b></DT>
+<DD>Verbose module output</DD>
+<DT><b>--quiet</b></DT>
+<DD>Quiet module output</DD>
+</DL>
+
+<h3>Parameters:</h3>
+<DL>
+<DT><b>DTM</b>=<em>string</em></DT>
+<DD>Input DTM raster map</DD>
+
+<DT><b>W_S_PROFILE</b>=<em>string</em></DT>
+<DD>Input ASCII file of the water surface profile</DD>
+
+<DT><b>RIVER</b>=<em>string</em></DT>
+<DD>Input vector line map of river-axis</DD>
+
+<DT><b>FLOODING_MAP</b>=<em>string</em></DT>
+<DD>Output: name of flooding map</DD>
+
+<DT><b>DOUBT_MAP</b>=<em>string</em></DT>
+<DD>Output: name of doubful surface areas</DD>
+
+<DT><b>PROFILE_T100</b>=<em>string</em></DT>
+<DD>Input ASCII file with water-depht for return period T > 100 years</DD>
+
+<DT><b>(null)</b>=<em>float</em></DT>
+<DD>Input delta_y to find the boundaries of the main channel [default value 0.5 m]</DD>
+
+<DT><b>delta_x</b>=<em>float</em></DT>
+<DD>Input delta_x to find the boundaries of the main channel [default value 3.5 m]</DD>
+
+<DT><b>river_boundary</b>=<em>string</em></DT>
+<DD>Output vector boundaries of the main channel</DD>
+
+<DT><b>boundary_type</b>=<em>string</em></DT>
+<DD>Format of vector boundaries of the main channel</DD>
+<DD>Options: <em>line,points</em></DD>
+<DD>Default: <em>points</em></DD>
+
+<DT><b>res_B</b>=<em>integer</em></DT>
+<DD>Input value: resolution B [default value 10 m]</DD>
+
+<DT><b>res_C</b>=<em>integer</em></DT>
+<DD>Input value: resolution C [default value 20 m]</DD>
+
+<DT><b>MAP_W_S_PROFILE</b>=<em>string</em></DT>
+<DD>Output vector point map of the water surface profile</DD>
+
+</DL>
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html>
+<head>
+</head>
+<body>
+<h2>DESCRIPTION</h2>
+
+
+<em><b>r.inund.fluv</b></em> - This command allows to obtain a fluvial potentially inundation map given a high-resolution DTM of the area surrounding the river and a water surface profile calculated through an 1-D hydrodinamic model.<br>
+
+The implemented procedure has innovating characteristics; even if it remains substantially one-dimensional, it takes into account the two-dimensionality of territory and inundation phenomena, adducing hypotheses that let to correct many errors typical of one-dimensional usually employed procedure. With respect to a two-dimensional model, it has advantage that it needs a lower computational effort, that allows to apply it to river reaches very long (of the order of 100km).<br>
+
+The single phases in which the procedure is divided, are here summarized, even if the whole procedure is excecuted automatically. <br>
+
+- In the first phase, the value of water level in the nearest point of fluvial axis is assigned to each pixel of terrain, through the creation of Thiessen polygons (<EM><A HREF="r.surf.idw.html">r.surf.idw</A></EM> setting <em>npoints</em>=1) <em>(default resolution value = 10 m)</em>. Then the procedure makes a comparison between elevation of each pixel and water level, and defines pixels characterized by elevation lower than water level, at hazard. <br>
+- In the second phase, the procedure removes all the areas previously defined at hazard but not connected with the river axis, i.e. surrounded by terrain not at hazard, neglecting infiltration or underground rivers (<EM><A HREF="v.select.html">v.select</A></EM>). <br>
+- In the third phase, the hypotesis is that water diffuses from river to the surrounding areas only in direction perpendicular to the river axis. Through an implemented fortran code (<em>clean_inundation.f90</em>), the procedure individuates the pixels, considered at hazard at the end of phase 2, inundated for sure by water; hence it dries the ones, considered at hazard at the end of phase 2, not reached by water because protected by levees or small hills posed along the perpedicular path between the pixel and the river axis <em>(default resolution value = 20 m)</em>. <br>
+- In the fourth phase, the hypothesis is that water, outside the main channel, moves along the maximum terrain slope direction. First, through an implemented fortran code (<em>find_main_channel.f90</em>), the procedure individuates the "boundaries" of the main channel. Then, through an implemented fortran code (<em>2d_path.f90</em>), the procedure individuates the water path outside the main channel along the maximum terrain slope direction <em>(default resolution value: 20 m)</em>. Finally, through an implemented fortran code (<em>correction_from_path.f90</em>), the procedure individuates the pixels, dried at the end of phase 3, connected to a water path, and defines them at hazard <em>(default resolution value: 10 m)</em>. <br>
+- In the fifth and last phase, the final potentially inundated map is defined as sum of areas at hazard in the third and fourth phases.
+<br>
+<br>
+<em><b><font size = "5">Use</b></em></font size = "5">
+<br>
+<p>
+<em><b>Requested input and output:</b></em> <br>
+The command requires some input:<br>
+
+1 - a Digital Terrain Model (DTM, DEM or DSM) of the area surrounding the river; an high-resolution DTM (e.g. 5 x 5 m or higher) is required to describe levees and other human infrastructures on the floodplain. <br>
+2 - an ASCII file describing the water surface profile along the channel axis, in this format (e.g.):
+<div class="code"><PRE>
+
+411815.62874469644 4944870.642304279 197. 104
+411848.8162966241 4944958.868788462 196.96 103.933*
+411882.0010118869 4945047.096337882 196.89 103.866*
+411915.188102181 4945135.322989558 196.82 103.8*
+ ...
+</PRE></div>
+where the first column is East-coordinate, the second is North-coordinate, the third is the water level in each cross-section, and the last one is the name of the cross-section. The distance between two cross-sections is suggested to be less than 100 meters. <br>
+3 - a vector line describing the river axis from upstream to dowstream. <br>
+<br>
+The output are: <br>
+4 - the fluvial potentially inundation map; <br>
+5 - the area in doubt, i.e. areas which could be at flooding hazard because of their elevation, but aren't defined at hazard at the end of the procedure.
+<br>
+<p>
+<em><b><font color="red">WARNING </b></em></font> -
+The original resolution of DTM must be high. If the pixel dimension is greater than 5 m, the command informs the user with a <em>warning message</em>; however the procedure runs.
+<br><br>
+ <em><b>Optional input and output:</b></em> <br>
+The command has some input parameters that the user can choose in the <em>Optional parameters & maps </em> interface. They refer to intermediate phases of the procedure, described in detail in reference.<br>
+1b - It's possible to set another water surface profile associated with a return period of 100 years or greater, from which the procedure starts to find automatically the "boundaries" of main channel. The file format is as mentioned above for the input file 2. <br>
+2b - It's possible to set a <em>delta_y</em> different from the default value (0.5 m), as the step with which the procedure increases the water level inside the main channel, to find automatically the "boundaries" of main channel itself.
+3b - It's possible to set a <em>delta_x</em> different from the default value (1.5*DTM_resolution in meters), i.e. the threshold value of the distance between the points individuated on a bank in two consecutive steps, to exceed so that the procedure may define the "boundary" of main channel in a pixel between such two points. <br>
+4b - The user may view the "boundaries" of main channel, to verify the working of the procedure, especially if <em>delta_y</em> and <em>delta_x</em> were setted different from the default values. <br>
+5b - 6b - The procedure performs some calculations using resolution different from the DTM one. The default values of 10 and 20 meters, used in different steps of the procedure, are a compromise between computational speed and consistency of output map. It's suggested that only the expert user modifies such values.
+<br>
+7b - The user can view a vector point map of the water surface profile along the river axis.
+<br>
+<br>
+<em>(GRASS Shell Script)</em>
+
+
+<h2>AUTHORS</h2>
+Roberto Marzocchi, University of Genoa, Italy <a href="mailto:roberto.marzocchi at gmail.com"> e-mail</a><br>
+Bianca Federici, University of Genoa, Italy <a href="mailto:bianca.federici at unige.it"> e-mail </a><br>
+Domenico Sguerso, University of Genoa, Italy <a href="mailto:domenico.sguerso at unige.it"> e-mail </a><br><p>
+
+<h2>REFERENCES</h2>
+- Federici B. & Sguerso D. (2008). Procedura automatica per la creazione di mappe di potenziale inondazione fluviale. In fase di revisione per la pubblicazione sul Bollettino SIFET. <br>
+- Pdf presentation of the work at the "IX Meeting degli Utenti Italiani di GRASS - GFOSS": <a href="http://www.grassmeeting2008.unipg.it/?q=node/9/"> web-page </a></em> <br>
+
+<p><i> Last changed: $21 april 2008 12:54:55 CET $</i></p>
+</body>
+</html>
Added: grass-addons/raster/r.inund.fluv/r.inund.fluv_62
===================================================================
--- grass-addons/raster/r.inund.fluv/r.inund.fluv_62 (rev 0)
+++ grass-addons/raster/r.inund.fluv/r.inund.fluv_62 2008-05-02 12:24:52 UTC (rev 31204)
@@ -0,0 +1,498 @@
+#!/bin/sh
+############################################################################
+#
+# MODULE: r.inund.fluv v.1 for GRASS 6.2
+# AUTHOR(S): Roberto Marzocchi, Bianca Federici, Domenico Sguerso
+# PURPOSE: Creates a fluvial inundation map given an high-resolution dtm and a water surface profile
+# COPYRIGHT: (C) 2008 by the GRASS Development Team
+#
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+#############################################################################
+
+#%Module
+#% description: Creates a fluvial inundation map given an high-resolution dtm and a water surface profile
+#% keywords: Automatic procedure to compute a fluvial inundation map
+#%End
+#%option
+#% key: DTM
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Input raster map (DTM or DEM)
+#% required : yes
+#%end
+#%option
+#% key: W_S_PROFILE
+#% type: string
+#% gisprompt: old_file,file,input
+#% description: Input ASCII file of the water surface profile
+#% required : yes
+#%end
+#%option
+#% key: RIVER
+#% type: string
+#% gisprompt: old,vector,vector
+#% description: Input vector map of river-axis
+#% required : yes
+#%end
+#%option
+#% key: FLOODING_MAP
+#% type: string
+#% gisprompt: new,cell,raster
+#% description: Output: name of flooding map
+#% required : yes
+#%end
+#%option
+#% key: DOUBT_MAP
+#% type: string
+#% gisprompt: new,cell,raster
+#% description: Output: name of doubful surface areas.
+#% required : no
+#%end
+#%option
+#% key: PROFILE_T100
+#% type: string
+#% gisprompt: old_file,file,input
+#% description: Input ASCII file with water-depht for return period T > 100 years
+#% required : no
+#%end
+#%option
+#% key: delta_y
+#% type: double
+#% gisprompt: new
+#% description: Input delta_y to find the boundaries of the main channel [default value 0.5 m]
+#% required : no
+#%end
+#%option
+#% key: delta_x
+#% type: double
+#% gisprompt: new
+#% description: Input delta_x to find the boundaries of the main channel [default value DTM_resolution * 1.5 m ]
+#% required : no
+#%end
+#%option
+#% key: river_boundary
+#% type: string
+#% gisprompt: new,vector,vector
+#% description: Output vector boundaries of the main channel
+#% required : no
+#%end
+#%option
+#% key: boundary_type
+#% type: string
+#% description: Format of vector boundaries of the main channel
+#% options: line,points
+#% answer: points
+#% required: no
+#%end
+#%option
+#% key: res_B
+#% type: integer
+#% gisprompt: new,cell,raster
+#% description: Input value: resolution B [default value 10 m]
+#% required : no
+#%end
+#%option
+#% key: res_C
+#% type: integer
+#% gisprompt: new,cell,raster
+#% description: Input value: resolution C [default value 20 m]
+#% required : no
+#%end
+#%option
+#% key: MAP_W_S_PROFILE
+#% type: string
+#% gisprompt: new,vector,vector
+#% description: Output vector map of the water surface profile
+#% required : no
+#%end
+
+if [ -z "$GISBASE" ]
+then
+ echo ""
+ echo "You must be in GRASS GIS to run this program"
+ echo ""
+ exit 1
+fi
+
+if [ "$1" != "@ARGS_PARSED@" ] ; then
+ exec g.parser "$0" "$@"
+fi
+
+$GRASS_VERBOSE
+
+#request control & test (only for imput map)
+if [ -z "$GIS_OPT_DTM" ] ; then
+ echo "ERROR: DTM not specified"
+ exit 1
+fi
+#g.findfile element=cell file="$GIS_OPT_DTM" > /dev/null # test if a map exists
+#if [ $? -eq 0 ]; then
+# echo "Imput map not found" 1>&2
+# exit 1
+#fi
+g.region save=dtm_res2 --overwrite
+g.region rast="$GIS_OPT_DTM" save=verifica
+g.region -p > .resolution.txt
+g.region region=dtm_res2
+exec 6<&0
+exec < .resolution.txt
+read a1 b1
+read a2 b2
+read a3 b3
+read a4 b4
+read a5 b5
+read a6 b6
+read a7 b7
+read a8 b8
+read a9 b9
+exec 0<&6 6<&-
+res=$b9
+if [ "$b9" -gt 5 ]
+then
+ echo "WARNING - the DTM resolution it'snt good for this application"
+else
+ echo "The DTM resolution is $res x $res meters"
+fi
+
+
+rm .resolution.txt
+
+
+if [ -z "$GIS_OPT_W_S_PROFILE" ] ; then
+ echo "ERROR: file with depht not specified"
+ exit 1
+fi
+
+
+
+if [ -z "$GIS_OPT_RIVER" ] ; then
+ echo "ERROR: vector with river-axis not specified"
+ exit 1
+fi
+
+
+if [ -z "$GIS_OPT_FLOODING_MAP" ] ; then
+ echo "ERROR: name of flooding map not specified"
+ exit 1
+fi
+
+
+#optional control: assignement default value
+if [ -n "$GIS_OPT_res_B" ]
+then
+ res10="$GIS_OPT_res_B"
+else
+ res10=10
+fi
+
+if [ -n "$GIS_OPT_res_C" ]
+then
+ res20="$GIS_OPT_res_C"
+else
+ res20=20
+fi
+
+if [ -n "$GIS_PROFILE_T100" ]
+then
+ quote100="$GIS_OPT_PROFILE_T100"
+else
+ quote100="$GIS_OPT_W_S_PROFILE"
+fi
+
+if [ -n "$GIS_OPT_delta_x" ]
+then
+ deltax="$GIS_OPT_delta_x"
+else
+ deltax=$(echo "$res * 1.5" | bc)
+fi
+
+if [ -n "$GIS_OPT_delta_y" ]
+then
+ deltay="$GIS_OPT_delta_y"
+else
+ deltay="0.5"
+fi
+
+g.remove rast=MASK
+cartella_nascosta=".temp_grass_script"
+# do not use any space in name of cartella_nascosta
+
+
+rm -r -f ~/$cartella_nascosta
+mkdir ~/$cartella_nascosta
+cp "$GISBASE/fortran_code"/find_main_channel.exe ~/$cartella_nascosta
+cp "$GISBASE/fortran_code"/clean_inundation.exe ~/$cartella_nascosta
+cp "$GISBASE/fortran_code"/2d_path.exe ~/$cartella_nascosta
+cp "$GISBASE/fortran_code"/correction_from_path.exe ~/$cartella_nascosta
+
+#variables
+dtm="$GIS_OPT_DTM"
+cp "$GIS_OPT_W_S_PROFILE" ~/$cartella_nascosta/profilo
+v.in.ascii -z input="$GIS_OPT_W_S_PROFILE" output=profilo format=point fs=" " skip=0 x=1 y=2 z=3 cat=0 --overwrite
+centerline="$GIS_OPT_RIVER"
+inondazione="$GIS_OPT_FLOODING_MAP"
+
+
+#calculation of the "boundaries" of the main channel
+g.region nsres=$res ewres=$res save=dtm_res2 --overwrite
+g.region -p >> ~/$cartella_nascosta/"region2.txt"
+r.out.ascii -h input=$dtm output=$cartella_nascosta/"dtm2x2" null=0
+
+#fortran code run
+cd
+cp "$quote100" ~/$cartella_nascosta/quote_100
+cd $cartella_nascosta
+echo $res $deltax $deltay > parameter.txt
+echo "'"region2.txt"'" "'"dtm2x2"'" "'"limiti_alveo_vect"'" "'"quote_100"'" > nomefile.txt
+echo "Fortran code that automatically find bed river"
+./find_main_channel.exe
+rm nomefile.txt
+rm parameter.txt
+#TEST FORTRAN CODE RUN
+if [ -z "limiti_alveo_vect" ]
+then
+ echo "ERROR IN FORTRAN CODE - see fortran error message and try understand the problem"
+ exit 1
+fi
+cd
+v.in.ascii --o input=$cartella_nascosta/limiti_alveo_vect output=limiti_alveo format=point fs=' ' skip=0 x=1 y=2 cat=0
+
+
+#export of 20x20 dtm, region20.txt and region10.txt
+g.region nsres=$res20 ewres=$res20 save=dtm_res20 --overwrite
+g.region -p >> ~/$cartella_nascosta/region20.txt
+r.out.ascii -h input=$dtm output=$cartella_nascosta/dtm20x20 null=0
+g.region nsres=$res10 ewres=$res10 save=dtm_res10 --overwrite
+g.region -p >> ~/$cartella_nascosta/region10.txt
+
+
+echo "STEP 1"
+###################################################################################################################
+# STEP 1 - Thiessen interpolation
+###################################################################################################################
+# don't use low resolution [default 10m]
+r.mask input=$dtm maskcats=*
+g.region region=dtm_res10
+v.surf.idw input=profilo output=quote_punti_adiacenti npoints=1 layer=1 column=dbl_3
+g.region region=dtm_res2
+r.mapcalculator amap=quote_punti_adiacenti bmap=$dtm formula="if(A>=B)" outfile=inondazione_step1 help=-
+
+
+echo "STEP 2"
+###################################################################################################################
+# STEP 2 - remove "lakes"
+###################################################################################################################
+r.to.vect -s input=inondazione_step1 output=inondazione_step1 feature=area --overwrite
+v.select ainput=inondazione_step1 atype=line,area alayer=1 binput=$centerline btype=line blayer=1 output=inondazione_step2 operator=overlap --overwrite
+v.to.rast input=inondazione_step2 output=inondazione_step2 use=val layer=1 value=1 rows=4096 --overwrite
+
+
+echo "STEP 3"
+###################################################################################################################
+# STEP 3 - looking at the minimum path from wet pixel to the channel axis
+# if the dtm elevation is higher than the water elevation in the pixel ==> pixel becomes dry
+# if the dtm elevation is always lower than the water elevation in the pixel ==> pixel remains wet
+###################################################################################################################
+# export data
+#20x20
+g.region region=dtm_res20
+r.out.ascii -h input=inondazione_step2 output=$cartella_nascosta/inondazione_step2 null=0
+r.out.ascii -h input=quote_punti_adiacenti output=$cartella_nascosta/quote_punti_adiacenti dp=2 null=0
+
+# FORTRAN code
+echo "CODE FORTRAN OF STEP 3"
+cd
+cd $cartella_nascosta
+echo "'"inondazione_step2"'" "'"profilo"'" "'"inondazione_nuova.txt"'" "'"quote_punti_adiacenti"'" "'"region20.txt"'" "'"region2.txt"'" "'"dtm2x2"'" > nomefile.txt
+./clean_inundation.exe
+rm nomefile.txt
+
+#TEST FORTRAN CODE RUN
+if [ -z "inondazione_nuova.txt" ]
+then
+ echo "ERROR IN FORTRAN CODE - see fortran error message and try understand the problem"
+ exit 1
+fi
+cd
+
+# output file
+r.in.ascii input=$cartella_nascosta/inondazione_nuova.txt output=inondazione_step3 'mult=1.0 or read from header' nv=0 --overwrite
+
+
+echo "STEP 4"
+###################################################################################################################
+# STEP 4 - calcolation of the 2d_path outside the main channel
+###################################################################################################################
+r.out.ascii -h input=inondazione_step3 output=$cartella_nascosta/inondazione_step3 null=0
+
+echo "2D MODEL OF DIFFUSION OF FLOODING INONDATION OUTSIDE THE MAIN CHANEL"
+cd
+cd $cartella_nascosta
+echo "'"profilo"'" "'"inondazione_step3"'" "'"reticolo"'" "'"region20.txt"'" "'"dtm20x20"'" "'"limiti_alveo_vect"'" "'"region2.txt"'" "'"dtm2x2"'" "'"bordi_alveo"'" > nomefile.txt
+./2d_path.exe
+rm nomefile.txt
+
+#TEST FORTRAN CODE RUN
+if [ -z "reticolo" ]
+then
+ echo "ERROR IN FORTRAN CODE - see fortran error message and try understand the problem"
+ exit 1
+fi
+cd
+
+
+#r.in.ascii input=$cartella_nascosta/bordi_alveo output=bordi_alveo 'mult=1.0 or read from header' nv=0 --overwrite
+#r.thin input=bordi_alveo output=bordi_alveo_thin iterations=200 --overwrite
+#r.to.vect input=bordi_alveo_thin output=bordi_alveo feature=line --overwrite
+v.in.ascii -n input=$cartella_nascosta/bordi_alveo output=bordi_alveo format=standard fs=' ' skip=0 x=1 y=2 z=0 cat=0
+r.in.ascii input=$cartella_nascosta/reticolo output=reticolo 'mult=1.0 or read from header' nv=* --overwrite
+r.thin input=reticolo output=reticolo_thin iterations=200 --overwrite
+r.to.vect input=reticolo_thin output=reticolo feature=line --overwrite
+g.remove rast=reticolo_thin
+
+
+# 10x10 resolution
+---> 1 pò più lento ma più preciso!
+g.region region=dtm_res10
+
+# difference between maps of step1 and step2
+r.mapcalculator amap=inondazione_step3 formula='isnull(A)' outfile=inondazione_step3_isnull help=-
+r.mapcalculator amap=inondazione_step3_isnull formula='if(A,0,1)' outfile=inondazione_step3 help=-
+r.mapcalculator amap=inondazione_step2 bmap=inondazione_step3 formula='A-B' outfile=diff_inond help=-
+r.mapcalculator amap=diff_inond formula='if(A,1,0,0)' outfile=diff_inond help=- #remove the -1 value resulted from A-B subtraction
+r.null map=diff_inond setnull=0
+r.to.vect input=diff_inond output=diff_inond1 feature=area --overwrite
+#to use the v.select command the line must have a layer associated ==> export of vector map in dxf format and then import the dxf file
+v.out.dxf input=reticolo output=$cartella_nascosta/reticolo.dxf
+v.in.dxf -1 input=$cartella_nascosta/reticolo.dxf output=reticolo --overwrite
+
+#only 2d_path overlapped the surfaces to look
+v.select ainput=reticolo atype=line alayer=1 binput=diff_inond1 btype=area blayer=1 output=reticolo_interno operator=overlap --overwrite
+
+# creation of points by 2d path outside the main channel
+v.to.points input=reticolo_interno type=point,line,boundary,centroid output=punti_reticolo_interno llayer=1 dmax=30 --overwrite
+
+# convert 2d points to 3d points with the elevation of the water surface profile at the starting point of the path in the channel axis (or with the elevation of the nearest water surface profile)
+g.region region=dtm_res20
+v.drape input=punti_reticolo_interno type=point,centroid,line,boundary,face,kernel rast=reticolo method=nearest output=punti3d_reticolo --overwrite
+#v.drape input=punti_reticolo_interno type=point,centroid,line,boundary,face,kernel rast=quote_punti_adiacenti method=nearest output=punti3d_reticolo --overwrite
+
+
+g.region region=dtm_res10
+
+#export of 3d points
+v.out.ascii input=punti3d_reticolo output=$cartella_nascosta/punti3d_reticolo format=point dp=2
+r.out.ascii -h input=diff_inond output=$cartella_nascosta/diff_inond null=0
+
+
+
+echo "FORTRAN CODE OF STEP 4"
+cd
+cd $cartella_nascosta
+echo "'"diff_inond"'" "'"punti3d_reticolo"'" "'"correzione"'" "'"region10.txt"'" "'"dtm2x2"'" "'"region2.txt"'" > nomefile.txt
+./correction_from_path.exe
+rm nomefile.txt
+
+#TEST FORTRAN CODE RUN
+if [ -z "correzione" ]
+then
+ echo "ERROR IN FORTRAN CODE - see fortran error message and try understand the problem"
+ exit 1
+fi
+cd
+
+
+r.in.ascii input=$cartella_nascosta/correzione output=inondazione_step5 'mult=1.0 or read from header' nv=0 --overwrite
+
+# replace null cells with 0 cells
+g.region region=dtm_res2
+r.mapcalculator amap=inondazione_step5 formula='isnull(A)' outfile=inondazione_step5_isnull help=-
+r.mapcalculator amap=inondazione_step5_isnull formula='if(A,0,1)' outfile=inondazione_step5 help=-
+r.mapcalculator amap=inondazione_step5 bmap=inondazione_step3 formula='A+B' outfile=inondazione help=-
+
+
+echo "STEP 5"
+###################################################################################################################
+# STEP 5 - EXPORT OUTPUT FILES
+###################################################################################################################
+g.rename rast=inondazione,inondazione1 --overwrite
+r.null map=inondazione1 setnull=0 --overwrite
+r.to.vect -s input=inondazione1 output=inondazione1 feature=area --overwrite
+v.patch input=centerline,punti3d_reticolo output=overlap --overwrite
+v.select ainput=inondazione1 atype=area alayer=1 binput=overlap btype=line blayer=1 output=inondazione operator=overlap --overwrite
+v.to.rast input=inondazione output=$inondazione use=val layer=1 value=1 rows=4096 --overwrite
+
+#COLOURS OF INUNDATION RASTER MAP
+rm -f ".colore_inondazione.file"
+echo '1 blue' >> ".colore_inondazione.file"
+echo '0 grey' >> ".colore_inondazione.file"
+echo 'nv white' >> ".colore_inondazione.file"
+echo 'end' >> ".colore_inondazione.file"
+cat ".colore_inondazione.file" | r.colors map=$inondazione color=rules
+rm ".colore_inondazione.file"
+
+
+if [ -n "$GIS_OPT_DOUBT_MAP" ]
+then
+ # It's possible to create a map with doubful surface areas
+ dubbio="$GIS_OPT_DOUBT_MAP"
+ r.mapcalculator amap=$inondazione formula='isnull(A)' outfile=temp help=-
+ r.mapcalculator amap=temp formula='if(A,0,1)' outfile=temp help=-
+ r.mapcalculator amap=inondazione_step2 bmap=temp formula='A-B' outfile=$dubbio help=-
+ r.null map=$dubbio setnull=0
+ #COLOURS OF DOUBFUL RASTER MAP
+ rm -f .colore_dubbio.file
+ echo '1 red' >> .colore_dubbio.file
+ echo '0 grey' >> .colore_dubbio.file
+ echo 'nv white' >> .colore_dubbio.file
+ echo 'end' >> .colore_dubbio.file
+ cat .colore_dubbio.file | r.colors map=$dubbio color=rules
+ rm .colore_dubbio.file
+ g.remove rast=temp
+fi
+
+
+g.remove rast=MASK
+g.remove rast=inondazione_step1,inondazione_step2,inondazione_step3,reticolo,quote_punti_adiacenti
+g.remove rast=inondazione_step3_isnull,diff_inond,inondazione_step5,inondazione_step5_isnull,inondazione1
+g.remove vect=inondazione_step1,inondazione_step2,reticolo
+g.remove vect=diff_inond1,reticolo_interno,inondazione1,overlap
+
+if [ "$inondazione" == "inondazione" ]
+then
+ echo " "
+else
+ g.remove vect=inondazione
+fi
+
+if [ -n "$GIS_OPT_MAP_W_S_PROFILE" ]
+then
+ g.rename vect=profilo,"$GIS_OPT_MAP_W_S_PROFILE"
+else
+ g.remove vect=profilo
+fi
+
+
+if [ -n "$GIS_OPT_river_boundary" ]
+then
+ if [ "$GIS_OPT_boundary_type" = "line" ]
+ then
+ g.rename vect=bordi_alveo,"$GIS_OPT_river_boundary"
+ g.remove vect=limiti_alveo
+ else
+ g.rename vect=limiti_alveo,"$GIS_OPT_river_boundary"
+ g.remove vect=bordi_alveo
+ fi
+else
+ g.remove vect=bordi_alveo,limiti_alveo
+fi
+
+rm -r ~/$cartella_nascosta
+
+
+
+exit 0
+
More information about the grass-commit
mailing list