C Program to set up pointers from DEM elevation matrix file and adjust C elevations of pits by increasing their elevation till they drain, C i.e. pointers can be assigned consistently. C C CREATED BY DAVID G TARBOTON C C Input is read from file BRANCH.IN which has the following structure. C Record 1: File name of raw elevation data file. C Record 2: File name of pointer file to be output. C Record 3: File name of area file (Not used). C Record 4: File name of adjusted elevation data file to be output. C Subsequent records: Number of subdivisions in X and Y directions for C matrix to be divided into in first sweeps to speed up adjustment of C elevations. C MEANING OF POINTERS IS ------------- C I 4 I 3 I 2 I C 0 = POINTS TO SELF ------------- C I.E. UNRESOLVED I 5 I 0 I 1 I C -1= BOUNDARY PIXEL ------------- C I 6 I 7 I 8 I C ------------- C C C This program uses subroutines from DEMUTIL.FOR, the set of subroutines C for I/O of binary matrix files. C C Program main PARAMETER(IGRIDY=2404,IGRIDX=2404) PARAMETER(IGY=IGRIDY/2+2,IGX=IGRIDX/2+2) INTEGER*2 ELEV(IGRIDY,IGRIDX),DIR(IGRIDY,IGRIDX), + ELEVP(IGY,IGX) integer*4 nx, ny, npx, npy, nxp, nyp real*4 dx, dy CHARACTER*80 DEMFILE,POINTFILE,AREAFILE,NEWFILE C C READ INPUT C OPEN(UNIT=11,form='formatted',FILE='branch.in',STATUS='OLD') READ(11,22) DEMFILE,POINTFILE,AREAFILE,NEWFILE 22 FORMAT(A80/A80/A80/A80) CALL demread (ELEV, DEMFILE, NX, NY, IGRIDY, DX, DY) if(nx.gt.igridx.or.ny.gt.igridy)then write(6,38) nx, igridx, ny, igridy 38 format(1x,'array dimensions too small'/ + 1x,'x',2i8,' y',2i8) stop endif C C LOOP HERE FOR NESTED PARTITIONS C 1 READ(11,*,END=99)NPX,NPY DO 2 IP=1,NPY DO 2 JP=1,NPX I1=MAX(((IP-1)*NY)/NPY,1) I2=(IP*NY)/NPY J1=MAX(((JP-1)*NX)/NPX,1) J2=(JP*NX)/NPX NXP=J2-J1+1 NYP=I2-I1+1 IF(NXP.LE.IGX.AND.NYP.LE.IGY)THEN C C LOAD ELEVATIONS INTO PARTITION ARRAY C DO 3 I=I1,I2 DO 3 J=J1,J2 ELEVP(I-I1+1,J-J1+1)=ELEV(I,J) 3 CONTINUE WRITE(6,100)IP,JP 100 FORMAT(1X,'PARTITION',2I5) CALL SETDF(ELEVP,DIR,IGX,IGY,NXP,NYP,DX,DY) C C WRITE ELEVATIONS BACK INTO MAIN ARRAY C DO 4 I=I1,I2 DO 4 J=J1,J2 ELEV(I,J)=ELEVP(I-I1+1,J-J1+1) 4 CONTINUE ELSE WRITE(6,*)'PARTITION TOO BIG',IP,JP,NXP,NYP ENDIF 2 CONTINUE GO TO 1 C C CALL SETDF FOR WHOLE AREA TO SMOOTH OFF JOINS C 99 CALL SETDF(ELEV,DIR,IGRIDX,IGRIDY,NX,NY,DX,DY) CALL demwrite (ELEV,NEWFILE,NX,NY,IGRIDY,DX,DY) CALL demwrite (DIR,POINTFILE,NX,NY,IGRIDY,DX,DY) END * * Subroutine to do the bulk of the work * SUBROUTINE SETDF(ELEV, DIR, IGRIDX, IGRIDY, NX, NY, DX, DY) INTEGER*2 ELEV(IGRIDY,IGRIDX),DIR(IGRIDY,IGRIDX), + IS(1000000),JS(1000000) real*4 FACT(8), dx, dy integer*4 nx, ny INTEGER*4 D1(8),D2(8) DATA D1/0,-1,-1,-1,0,1,1,1/ DATA D2/1,1,0,-1,-1,-1,0,1/ C C MEANING OF POINTERS IS ------------- C I 4 I 3 I 2 I C 0 = POINTS TO SELF ------------- C I.E. UNRESOLVED I 5 I 0 I 1 I C -1= BOUNDARY PIXEL ------------- C I 6 I 7 I 8 I C ------------- C C INITIALISE BOUNDARY POINTERS IN MATRIX DIR C DO 2 I=1,NX DIR(1,I)=-1 DIR(NY,I)=-1 2 CONTINUE DO 3 I=1,NY DIR(I,1)=-1 DIR(I,NX)=-1 3 CONTINUE IUP=0 c c---initialise internal pointers (-ve elevation indicates outside domain) c do 31 i=2,ny-1 do 31 j=2,nx-1 if(elev(i,j).le.0)then dir(i,j)=-1 else dir(i,j)=0 endif 31 continue C C TEST ALL INTERNAL ELEVATIONS AND SET POINTERS C C WRITE(6,*)'FACTORS' DO 21 K=1,8 FACT(K)=1./SQRT(D1(K)*DY*D1(K)*DY+D2(K)*D2(K)*DX*DX) C WRITE(6,*)K,FACT(K) 21 CONTINUE WRITE(6,*)'PROBLEM PIXELS' WRITE(6,*)' FLATS UNRESOLVED' C C TEST ALL INTERNAL ELEVATIONS AND SET POINTERS C 1 DO 4 I=2,NY-1 DO 4 J=2,NX-1 if(elev(i,j).gT.0) then CALL SET(ELEV,DIR,I,J,IGRIDY,igridx,FACT) end if 4 CONTINUE C C---FIRST FIXING PASS C ELIMINATE FLATS C STORE UNRESOLVED PIXELS IN A STACK N=0 DO 7 I=2,NY-1 DO 7 J=2,NX-1 IF(DIR(I,J).EQ.0)THEN N=N+1 IS(N)=I JS(N)=J ENDIF IF(N.GE.1000000)THEN WRITE(6,*)'ARRAYS NOT BIG ENOUGH 1', n C WRITE(6,*)I,J C WRITE(6,*)ELEV(I-1,J-1),ELEV(I-1,J),ELEV(I-1,J+1) C WRITE(6,*)ELEV(I,J-1),ELEV(I,J),ELEV(I,J+1) C WRITE(6,*)ELEV(I+1,J-1),ELEV(I+1,J),ELEV(I+1,J+1) RETURN ENDIF 7 CONTINUE NFLAT=N C---CALL ROUTINE TO MAKE FLATS DRAIN TO NEIGHBOUR'S CALL VDN(ELEV,DIR,IS,JS,IGRIDY,N) C C ANY UNRESOLVED PIXELS HERE ARE POOLS SO RAISE THEM AND START C AGAIN C IUP=IUP+1 DO 51 II=1,N I=IS(II) J=JS(II) C---DETERMINE ELEVATION OF LOWEST NEIGHBOUR ILN=MIN(ELEV(I,J+1),ELEV(I-1,J+1),ELEV(I-1,J),ELEV(I-1,J-1), & ELEV(I,J-1),ELEV(I+1,J-1),ELEV(I+1,J),ELEV(I+1,J+1)) C---INCREMENT IS 1 OR DIFFERENCE BETWEEN LOWEST NEIGHBOUR ELEV(I,J)=ELEV(I,J)+MAX(1,ILN-ELEV(I,J)) 51 CONTINUE WRITE(6,*)NFLAT,N C---TEST IF COMPLETE IF(N.GE.1)GO TO 1 ! LOOP BACK IF SOMETHING CHANGED RETURN END * * ROUTINE TO DRAIN UNRESOLVED PIXELS TOWARDS NEIGHBOURS * SUBROUTINE VDN(ELEV,DIR,IS,JS,IGRIDY,N) INTEGER*2 IS(*),JS(*),ELEV(IGRIDY,*),DIR(IGRIDY,*),ED INTEGER*2 DN(1000000) INTEGER*2 D1(8),D2(8),K DATA D1/0,-1,-1,-1,0,1,1,1/ DATA D2/1,1,0,-1,-1,-1,0,1/ C---INITIALISE NEW DIRECTIONS TO 0 1 DO 4 IP=1,N 4 DN(IP)=0 DO 2 K=1,7,2 DO 2 IP=1,N ED=ELEV(IS(IP),JS(IP))-ELEV(IS(IP)+D1(K),JS(IP)+D2(K)) IF(ED.GE.0.AND.DIR(IS(IP)+D1(K),JS(IP)+D2(K)).NE.0.AND. & DN(IP).EQ.0)DN(IP)=K ! LOGICAL EQUIVALENTS TO COND M'S BELOW C I1=CVMGP(1,0,ED) C---I1 IS 1 IF DROP IS .GE.0 C I3=CVMGZ(0,1,DIR(IS(IP)+D1(K),JS(IP)+D2(K))) C---I3 IS 0 UNLESS NEIGHBOUR HAS DIRECTION SET C DN(IP)=CVMGZ(i1*i3*k,dn(ip),DN(IP)) C---DN IS NEW DIRECTION POINTER 2 CONTINUE DO 12 K=2,8,2 DO 12 IP=1,N ED=ELEV(IS(IP),JS(IP))-ELEV(IS(IP)+D1(K),JS(IP)+D2(K)) IF(ED.GE.0.AND.DIR(IS(IP)+D1(K),JS(IP)+D2(K)).NE.0.AND. & DN(IP).EQ.0)DN(IP)=K ! LOGICAL EQUIVALENTS TO COND M'S BELOW C I1=CVMGP(1,0,ED) C---I1 IS 1 IF DROP IS .GE.0 C I2=CVMGZ(K,DIR(IS(IP),JS(IP)),DIR(IS(IP),JS(IP))) C---I3 IS 0 UNLESS NEIGHBOUR HAS DIRECTION SET C DN(IP)=CVMGZ(i1*i3*k,dn(ip),DN(IP)) C---DN IS NEW DIRECTION POINTER 12 CONTINUE NI=0 DO 3 IP=1,N IF(DN(IP).GT.0)THEN DIR(IS(IP),JS(IP))=DN(IP) ELSE NI=NI+1 IS(NI)=IS(IP) JS(NI)=JS(IP) ENDIF 3 CONTINUE C WRITE(6,*)'IN VDN: NI, N ',NI,N IF(NI.LT.N)THEN N=NI GO TO 1 ENDIF N=NI RETURN END * * SUBROUTINE TO SET POINTERS IN DIRECTION OF STEEPEST DECENT * SUBROUTINE SET(ELEV,DIR,I,J,IGRIDY,igridx,FACT) INTEGER*2 ELEV(IGRIDY,1),DIR(IGRIDY,igridx) real*4 SLOPE(8),FACT(8) INTEGER*4 D1(8),D2(8) DATA D1/0,-1,-1,-1,0,1,1,1/ DATA D2/1,1,0,-1,-1,-1,0,1/ DO 2 K=1,8 SLOPE(K)=FACT(K)*FLOAT(ELEV(I,J)-ELEV(I+D1(K),J+D2(K))) 2 CONTINUE SMAX=0. DIR(I,J)=0 DO 3 K=1,8 IF(SLOPE(K).GT.SMAX)THEN SMAX=SLOPE(K) DIR(I,J)=K ENDIF 3 CONTINUE RETURN END