forrest@0: program realizer forrest@0: ! program MidPtFM2d forrest@0: ! forrest@0: ! fractal generator adapted from algorithm MidPointFM2D in the forrest@0: ! Science of Fractal Images, M.F. Barnsley, R. L. Devaney forrest@0: ! B. B. Mandelbrot, H.-O. Peitgen, D. Saupe, and R. F. Voss forrest@0: ! Springer-Verlag forrest@0: ! forrest@0: ! this program modified to combine multiple fractal topographies forrest@0: ! then normalize them relative to each other, forrest@0: ! and perform multiple horizontal slices forrest@0: ! to obtain requested p values for each habitat forrest@0: ! in a new synthetic realization of the landscape forrest@0: ! forrest@0: ! William W. Hargrove forrest@0: ! Paul M Schwartz forrest@0: ! 5/7/96 forrest@0: ! (423) 241-2748 forrest@0: ! forrest@0: use realizer_mod forrest@0: implicit none forrest@0: forrest@0: integer(i4) :: P(maxN*maxN,0:maxsp+3) forrest@0: real(r4) :: X(0:maxN,0:maxN) forrest@0: real(r4) :: prob(0:maxsp) forrest@0: real(r4) :: ssum forrest@0: real(r4) :: ztemp forrest@0: integer(i4) :: igreatest forrest@0: integer(i4) :: hslice(0:maxsp) forrest@0: real(r4) :: elev(maxN*maxN) forrest@0: real(r4) :: H(maxN,maxsp) forrest@0: integer(i4) :: II(0:maxN,0:maxN) forrest@0: integer(i4) :: kcnt(0:maxsp) forrest@0: integer(i4) :: ties(maxsp) forrest@0: integer(i4) :: ifcnt forrest@0: integer(i4) :: iend forrest@0: integer(i4) :: istart forrest@0: integer(i4) :: icnt forrest@0: real(r4) :: sigma forrest@0: real(r4) :: xmin, xmax, xmean, xss forrest@0: real(r4) :: range1 forrest@0: real(r4) :: zmean forrest@0: real(r4) :: pfract forrest@0: real(r4) :: xlow forrest@0: integer(i4) :: maxlevel forrest@0: integer(i4) :: numpainted forrest@0: integer(i4) :: numtied forrest@0: integer(i4) :: numblank forrest@0: integer(i4) :: ihighest forrest@0: integer(i4) :: ibiggest forrest@0: integer(i4) :: idummy forrest@0: integer(i4) :: idtype forrest@0: integer(i4) :: idumm forrest@0: integer(i4) :: isum, jsum forrest@0: integer(i4) :: ipixid forrest@0: integer(i4) :: kount forrest@0: integer(i4) :: i, j, k, l, m ! loop indices forrest@0: integer(i4) :: iv, iz forrest@0: integer(i4) :: N forrest@0: integer :: useed ! user seed value forrest@0: integer :: rsize ! length of array of seeds forrest@0: integer :: ierr ! error forrest@0: integer, allocatable :: iseed(:) ! array of random number seeds forrest@0: integer(i4) :: idispvec(0:maxN*maxN) forrest@0: character (len=1) :: ans forrest@0: character (len=28) :: maptitle forrest@0: character (len=60) :: mapfile forrest@0: character (len=10) :: mapcolor forrest@0: character (len=40) :: fmt ! For dynamic format strings in gfortran forrest@0: logical :: addition, wrap, gradient, slice, invert forrest@0: logical :: grassout, xpmout, probout, vizout forrest@0: logical :: constraint(0:maxsp) forrest@0: ! forrest@0: ! write execution script for next run forrest@0: open (9, file='input.scr',status='unknown') forrest@0: ! forrest@0: ! forrest@0: ! enter random number seed forrest@0: print 196 forrest@0: 196 format (t10,'Do you want to enter a random number seed ', & forrest@0: /t12,' or use the system clock? [clock]') forrest@0: read (5,97) ans forrest@0: 97 format (a1) forrest@0: if (ans .eq. 'Y' .or. ans .eq. 'y') then forrest@0: write (9,197) ans forrest@0: 197 format (a1,t20,"(enter seed or clock?)") forrest@0: print 100 forrest@0: 100 format (t10,'Enter random number seed ') forrest@0: read *, useed forrest@0: write (9,101) useed forrest@0: 101 format (i12,t20,"(random number seed)") forrest@0: else forrest@0: 21 useed = int(time8(),kind(useed)) forrest@0: if (useed .eq. 0) then forrest@0: print *, "useed was zero." forrest@0: go to 21 forrest@0: end if forrest@0: print 298 forrest@0: 298 format (t15,'Random number seed from system clock.'/) forrest@0: write (9,29) ans, useed forrest@0: 29 format (a1,t20,"(random seed from clock = ",i12,")") forrest@0: end if forrest@0: ! Set the random number seed forrest@0: call random_seed(size=rsize) forrest@0: allocate(iseed(rsize),stat=ierr) forrest@0: if (ierr /= 0) then forrest@0: print *,'Error allocating iseed(rsize)' forrest@0: stop forrest@0: end if forrest@0: iseed(1:rsize) = useed forrest@0: call random_seed(put=iseed) forrest@0: ! forrest@0: print 1017, useed forrest@0: 1017 format (t15,'Random number seed = ',i12/) forrest@0: ! forrest@0: ! forrest@0: print 512 forrest@0: 512 format (t10,'Do you want to visualize', & forrest@0: /t12, ' the synthetic landscape? ') forrest@0: read (5,522) ans forrest@0: 522 format (a1) forrest@0: write (9,57) ans forrest@0: 57 format (a1,t20,"(visualize the landscape?)") forrest@0: if (ans .eq. 'y' .or. ans .eq. 'Y') then forrest@0: vizout = .true. forrest@0: else forrest@0: vizout = .false. forrest@0: end if forrest@0: ! forrest@0: 1 print 225 forrest@0: 225 format (t10,'How many levels? (map is 2**levels on a side)') forrest@0: read *, maxlevel forrest@0: ! forrest@0: N = 2**maxlevel forrest@0: irowcols = N forrest@0: numcells = (N+1)*(N+1) forrest@0: if (N .gt. maxN) goto 1 forrest@0: write (9,102) maxlevel forrest@0: 102 format (i3,t20,"(Number of levels)") forrest@0: ! forrest@0: print 114 forrest@0: 114 format (t10,'Highest data category (excluding no-data)', & forrest@0: /t12,' in the realized synthetic landscape?') forrest@0: read *, ihabs forrest@0: write (9,199) ihabs forrest@0: 199 format (i3,t20,"(highest data category)") forrest@0: ! forrest@0: print 232 forrest@0: 232 format (t10,'Will you supply external probability maps', & forrest@0: /t12,' for any of these categories?') forrest@0: read (5,111) ans forrest@0: 111 format (a1) forrest@0: write (9,161) ans forrest@0: 161 format (a1,t20,"(supply external prob maps?)") forrest@0: ! forrest@0: do j = 0, ihabs forrest@0: constraint(j) = .false. forrest@0: end do forrest@0: ! forrest@0: isupply = 0 forrest@0: if (ans .eq. 'Y' .or. ans .eq. 'y') then forrest@0: ! forrest@0: print 134 forrest@0: 134 format (t10,'For how many categories in the synthetic', & forrest@0: ' landscape',/t12,' will you supply external probability', & forrest@0: ' maps?') forrest@0: read *, isupply forrest@0: write (9,189) isupply forrest@0: 189 format (i3,t20,"(number of external prob maps)") forrest@0: ! forrest@0: do i = 1, isupply forrest@0: print 216 forrest@0: 216 format (t10,'Constrain which map category?') forrest@0: read *, k forrest@0: write (9,799) k forrest@0: 799 format (i3,t20,"(constrain which category?)") forrest@0: ! set constraint flag for this category forrest@0: constraint(k) = .true. forrest@0: ! forrest@0: print 220 forrest@0: 220 format (t10,'Enter name of input map file (integer*4, ', & forrest@0: /t12,' delimited by single spaces): ') forrest@0: read (5,230) mapfile forrest@0: 230 format (a60) forrest@0: print 2030, mapfile forrest@0: 2030 format (t15,'Input map name: ',a60/) forrest@0: write (9,231) mapfile forrest@0: 231 format (a60) forrest@0: ! forrest@0: print 144 forrest@0: 144 format (t10, 'Invert probabilities (i.e., DEM)? ') forrest@0: read (5,321) ans forrest@0: 321 format (a1) forrest@0: write (9,327) ans forrest@0: 327 format (a1,t20,"(invert probabilities?)") forrest@0: if (ans .eq. 'y' .or. ans .eq. 'Y') then forrest@0: invert = .true. forrest@0: else forrest@0: invert = .false. forrest@0: end if forrest@0: ! forrest@0: ! load constraint maps in proper column of P array forrest@0: ! invert probabilities if necessary, i.e., a DEM forrest@0: Call InputMap(mapfile,P,N,k,invert) forrest@0: end do forrest@0: ! forrest@0: end if forrest@0: ! forrest@0: ! forrest@0: do iz = 0, ihabs forrest@0: if (.not. constraint(iz)) then forrest@0: do i = 1, maxlevel forrest@0: 118 print 159, i, iz forrest@0: 159 format (t10,'Input value of H for level ',i3,', category ',i3) forrest@0: read *, H(i,iz) forrest@0: print *, H(i,iz) forrest@0: write (9,155) H(i,iz), i, iz forrest@0: 155 format (f4.2,t20,"(value of H for level ",i2,& forrest@0: ", category", i3,")") forrest@0: if (H(i,iz) .lt. 0.0 .or. H(i,iz) .gt. 1.0) goto 118 forrest@0: end do forrest@0: end if forrest@0: 664 end do forrest@0: ! forrest@0: ! forrest@0: ssum = 0.0 forrest@0: do i = 1, ihabs forrest@0: 115 print 116, i forrest@0: 116 format (t10,'Probability for habitat ',i3) forrest@0: read *, prob(i) forrest@0: if (prob(i) .lt. 0.0) goto 115 forrest@0: if (prob(i) .gt. 1.0) goto 115 forrest@0: ssum = ssum + prob(i) forrest@0: if (ssum .gt. 1.0000) then forrest@0: print *, 'Total probability exceeds 1.0 - Re-enter' forrest@0: ssum = ssum - prob(i) forrest@0: goto 115 forrest@0: end if forrest@0: write (9,151) prob(i), i forrest@0: 151 format (f6.5,t20,"(probability of habitat type", i2,")") forrest@0: end do forrest@0: ! forrest@0: ! find greatest non-constraint p to make background forrest@0: ztemp = 0._r4 forrest@0: igreatest = 0 forrest@0: do i = 0, ihabs forrest@0: if (.not. constraint(i)) then forrest@0: if (prob(i) .gt. ztemp) then forrest@0: ztemp = prob(i) forrest@0: igreatest = i forrest@0: !print *, 'at i of ',i ,'new greatest is ', igreatest forrest@0: end if forrest@0: end if forrest@0: 662 end do forrest@0: ! forrest@0: ! zero category by difference forrest@0: prob(0) = 1.0 - ssum forrest@0: if (prob(0) .lt. 0.00001) prob(0) = 0.0 forrest@0: print *, 'prob(0) is ', prob(0) forrest@0: ! forrest@0: ! print *, 'igreatest is ', igreatest forrest@0: if (prob(0) .gt. prob(igreatest)) igreatest = 0 forrest@0: print *, 'Cat ', igreatest,' will be calculated by absence' forrest@0: ! forrest@0: ! The Fractal Realizer assumes that the user-requested category forrest@0: ! with the greatest coverage on the landscape (the highest requested p) forrest@0: ! is the background or matrix category. In order to give the other forrest@0: ! categories the greatest freedom for fractal tuning and to minimize forrest@0: ! contention ties, this majority background color is computed by absence, forrest@0: ! and all other categories are reclassed. Furthermore, categories which forrest@0: ! are requested but for which the requested p = 0.0 are not actually forrest@0: ! computed. forrest@0: ! forrest@0: ! calculate reclassed categories forrest@0: ! ihabs = 0 forrest@0: ! do i = 0, ihabsmax forrest@0: ! if (i .eq. igreatest) then forrest@0: ! ireclass(0) = i forrest@0: ! print *, 'Cat ', i,' reassigned to background' forrest@0: ! goto 215 forrest@0: ! end if forrest@0: ! if (prob(i) .eq. 0.0) then forrest@0: ! print *, 'Cat ', i,' requested but not used - skipping' forrest@0: ! goto 215 forrest@0: ! end if forrest@0: ! ihabs = ihabs + 1 forrest@0: ! ireclass(ihabs) = i forrest@0: ! print *, 'cat ', i,' will now be cat ', ihabs forrest@0: ! forrest@0: !215 end do forrest@0: ! forrest@0: ! forrest@0: ! do i = 0, ihabs forrest@0: ! print 2002, ireclass(i), i, prob(ireclass(i)) forrest@0: !2002 format ('Original category ', i2, ' re-mapped to ', & forrest@0: ! i2, ' with requested probability of ', f6.4) forrest@0: ! end do forrest@0: ! forrest@0: print 212 forrest@0: 212 format (t10,'Generate GRASS r.in.ascii file? ') forrest@0: read (5,222) ans forrest@0: 222 format (a1) forrest@0: write (9,27) ans forrest@0: 27 format (a1,t20,"(generate GRASS input file?)") forrest@0: if (ans .eq. 'y' .or. ans .eq. 'Y') then forrest@0: grassout = .true. forrest@0: else forrest@0: grassout = .false. forrest@0: end if forrest@0: ! forrest@0: print 122 forrest@0: 122 format (t10,'Generate xpm file of synthetic landscape? ') forrest@0: read (5,123) ans forrest@0: 123 format (a1) forrest@0: write (9,124) ans forrest@0: 124 format (a1,t20,"(xpm file of realized landscape?)") forrest@0: if (ans .eq. 'y' .or. ans .eq. 'Y') then forrest@0: xpmout = .true. forrest@0: else forrest@0: xpmout = .false. forrest@0: end if forrest@0: ! forrest@0: print 322 forrest@0: 322 format (t10,'Generate xpm file of fractal probability' & forrest@0: ' landscapes? ') forrest@0: read (5,323) ans forrest@0: 323 format (a1) forrest@0: write (9,324) ans forrest@0: 324 format (a1,t20,"(xpm file of probability landscapes?)") forrest@0: if (ans .eq. 'y' .or. ans .eq. 'Y') then forrest@0: probout = .true. forrest@0: else forrest@0: probout = .false. forrest@0: end if forrest@0: ! forrest@0: ! ***************************************************************** forrest@0: ! fill P array with a separate fractal probability map forrest@0: ! for each desired habitat in the landscape forrest@0: do iz = 0, ihabs forrest@0: if (constraint(iz)) goto 22 forrest@0: ! forrest@0: ! don't bother if igreatest or no cells requested of this category forrest@0: if (prob(iz) .eq. 0.0 .or. iz .eq. igreatest) then forrest@0: do i = 1, numcells forrest@0: P(i,iz) = 0 forrest@0: end do forrest@0: goto 22 forrest@0: end if forrest@0: ! forrest@0: ! print 102 forrest@0: !102 format (t10,'Input the value for sigma') forrest@0: ! read *, sigma forrest@0: sigma = 1. forrest@0: print 103 forrest@0: 103 format (t10,'Input the threshold for detection (xlow) ') forrest@0: read (5,119) xlow forrest@0: 119 format (g12.6) forrest@0: write (9,125) xlow forrest@0: 125 format (f12.6,t20,"(threshold for detection (xlow))") forrest@0: if (xlow .eq. 0) xlow = -1e30 forrest@0: print *, xlow forrest@0: print 104 forrest@0: 104 format (t10,'Wrap? ') forrest@0: read (5,145) ans forrest@0: 145 format (a1) forrest@0: write (9,198) ans forrest@0: 198 format (a1,t20,"(wrap?)") forrest@0: if (ans .eq. 'y' .or. ans .eq. 'Y') then forrest@0: wrap = .true. forrest@0: else forrest@0: wrap = .false. forrest@0: end if forrest@0: ! forrest@0: ! no gradient possible with wrap forrest@0: if (.not. wrap) then forrest@0: print 106 forrest@0: 106 format (t10,'Gradient? ') forrest@0: read (5,107) ans forrest@0: 107 format (a1) forrest@0: write (9,148) ans forrest@0: 148 format (a1,t20,"(gradient?)") forrest@0: if (ans .eq. 'y' .or. ans .eq. 'Y') then forrest@0: gradient = .true. forrest@0: ! forrest@0: print 1850 forrest@0: 1850 format (t10,'Gradient: Initial values for four map corners') forrest@0: print 1851 forrest@0: 1851 format (t10,'Northeast? ') forrest@0: read *, northeast forrest@0: print 1852, northeast forrest@0: 1852 format (5x,g12.6) forrest@0: write (9,135) northeast forrest@0: 135 format (f12.6,t20,"(Gradient: northeast)") forrest@0: ! forrest@0: print 1853 forrest@0: 1853 format (t10,'Southeast? ') forrest@0: read *, southeast forrest@0: print 1852, southeast forrest@0: write (9,136) southeast forrest@0: 136 format (f12.6,t20,"(Gradient: southeast)") forrest@0: ! forrest@0: print 1854 forrest@0: 1854 format (t10,'Southwest? ') forrest@0: read *, southwest forrest@0: print 1852, southwest forrest@0: write (9,137) southwest forrest@0: 137 format (f12.6,t20,"(Gradient: southwest)") forrest@0: ! forrest@0: print 1855 forrest@0: 1855 format (t10,'Northwest? ') forrest@0: read *, northwest forrest@0: print 1852, northwest forrest@0: write (9,138) northwest forrest@0: 138 format (f12.6,t20,"(Gradient: northwest)") forrest@0: ! forrest@0: else forrest@0: gradient = .false. forrest@0: end if forrest@0: ! forrest@0: end if !wrap = .false. forrest@0: ! forrest@0: ! print 111 forrest@0: !111 format (t10,'Slice into discrete habitats? ') forrest@0: ! read (5,110) ans forrest@0: !110 format (a1) forrest@0: ! if (ans .eq. 'y' .or. ans .eq. 'Y') then forrest@0: ! slice = .true. forrest@0: ! else forrest@0: ! slice = .false. forrest@0: ! end if forrest@0: ! forrest@0: slice = .false. forrest@0: ! forrest@0: ! print 104 forrest@0: ! 104 format (t10,'Addition ?') forrest@0: ! read (5,105) ans forrest@0: ! 105 format (a1) forrest@0: ! if (ans .eq. 'y' .or. ans .eq. 'Y') then forrest@0: ! addition = .true. forrest@0: ! else forrest@0: ! addition = .false. forrest@0: ! end if forrest@0: addition = .false. forrest@0: ! forrest@0: ! call 2d multi-fractal generator forrest@0: print *,"Calling fractal generator for category ", iz forrest@0: call mpfm2d(X, maxlevel, sigma, H, addition, wrap, gradient, iz) forrest@0: ! forrest@0: print *, "fractal generated" forrest@0: xmin = 1e30 forrest@0: xmax = -1e30 forrest@0: xmean = 0.0 forrest@0: xss = 0.0 forrest@0: icnt = 0 forrest@0: ! forrest@0: ! obtain xmin and xmax for range normalization forrest@0: ! calculate initial raw statistics forrest@0: ! while you're at it forrest@0: do i = 0,N forrest@0: do j = 0,N forrest@0: if (X(i,j) .lt. xmin) xmin = X(i,j) forrest@0: if (X(i,j) .gt. xmax) xmax = X(i,j) forrest@0: ! xmean = xmean + X(i,j) forrest@0: ! xss = xss + X(i,j)*X(i,j) forrest@0: ! icnt = icnt+1 forrest@0: end do forrest@0: end do forrest@0: ! xmean = xmean / icnt forrest@0: ! print 200, xmean, xss, xmin, xmax forrest@0: !200 format (/t10,'Hfract Summary'/ & forrest@0: ! t10,'Mean = ',g12.6 & forrest@0: ! /t10,'Sum Sq. = ',g12.6 & forrest@0: ! /t10,'Minimum = ',g12.6 & forrest@0: ! /t10,'Maximum = ',g12.6) forrest@0: ! forrest@0: ! rescale and normalize fractal map between 0 and 32767, inclusive forrest@0: ! forrest@0: range1 = xmax - xmin forrest@0: do i = 0,N forrest@0: do j = 0,N forrest@0: iv = (i*(N+1))+j+1 forrest@0: P(iv,iz) = ((X(i,j) - xmin)/range1) * 32767 forrest@0: P(iv,ihabs+1) = i forrest@0: P(iv,ihabs+2) = j forrest@0: ! write(10,*) P(iv,iz) forrest@0: end do forrest@0: end do forrest@0: print*,'done normalizing now' forrest@0: ! forrest@0: !c write GRASS ascii map forrest@0: ! write(10,791), N forrest@0: !791 format ('north: ',i4) forrest@0: ! write(10,792) forrest@0: !792 format ('south: 0') forrest@0: ! write(10,793), N forrest@0: !793 format ('east: ',i4) forrest@0: ! write(10,794) forrest@0: !794 format ('west: 0') forrest@0: ! write(10,795), N+1 forrest@0: !795 format ('rows: ',i4) forrest@0: ! write(10,796), N+1 forrest@0: !796 format ('cols: ',i4) forrest@0: ! do i = 0,N forrest@0: ! write(10,*) (P((i*(N+1))+j+1,iz), j=0,N) forrest@0: ! end do forrest@0: ! forrest@0: ! end of iz loop over all landscape categories in synthetic map forrest@0: 22 continue forrest@0: end do forrest@0: ! ********************************************************************** forrest@0: ! forrest@0: ! ihabs fractal maps now installed in P matrix: forrest@0: ! forrest@0: ! ihabs cols + 3 forrest@0: ! |---------------| forrest@0: ! | | forrest@0: ! | | forrest@0: ! | empty cells | forrest@0: ! mtresolved =>| | forrest@0: ! | | forrest@0: ! mtsplit =>|---------------| forrest@0: ! | | forrest@0: ! | | forrest@0: ! | | forrest@0: ! | | forrest@0: ! | singly- | forrest@0: ! | painted | forrest@0: ! | | forrest@0: ! | | forrest@0: ! | | forrest@0: ! | | forrest@0: ! | | forrest@0: ! itiedsplit =>|---------------| forrest@0: ! | | forrest@0: ! | | forrest@0: ! | ties | forrest@0: ! m =>| | forrest@0: ! | | forrest@0: ! numcells =>|_______________| forrest@0: ! forrest@0: ! ihabs+1 column carries x-coord of cell forrest@0: ! ihabs+2 column carries y-coord of cell forrest@0: ! ihabs+3 column carries number of ties for this cell forrest@0: ! forrest@0: ! forrest@0: ! initialize II array with background value for re-use as output map forrest@0: II(:,:) = 32767 forrest@0: ! forrest@0: ! set paint thresholds forrest@0: print *,"Calculating paint thresholds" forrest@0: do i = 0,ihabs ! over all map categories forrest@0: ! forrest@0: !c don't bother if no cells requested of this cat forrest@0: if (prob(i) .ne. 0.0) then forrest@0: ! forrest@0: ! istart = 0 forrest@0: ! print *, "numcells is ", numcells forrest@0: ! print *, "istart is ", istart forrest@0: ! print *, "i is ", i forrest@0: ! Call HeapSort(P,numcells,istart,i) ! sort P matrix by ith column forrest@0: ! forrest@0: ! forrest@0: ! next category to be painted goes to single vector for sorting forrest@0: idispvec(0) = 0 forrest@0: do l = 1, numcells forrest@0: idispvec(l) = P(l,i) forrest@0: end do forrest@0: ! forrest@0: ! call special vector version of Heapsort forrest@0: istart = 1 forrest@0: Call HS2(idispvec,numcells+1,istart) forrest@0: ! forrest@0: !c diagnostic print forrest@0: ! do k=16636,16646 forrest@0: ! print *, i, k, idispvec(k) forrest@0: ! end do forrest@0: ! forrest@0: ! set paint threshold for this map category forrest@0: hslice(i) = idispvec(numcells - (int((prob(i) * numcells)+.499))) forrest@0: ! for absolute constraint maps forrest@0: if (hslice(i) .eq. 32767 .and. constraint(i)) hslice(i) = 32766 forrest@0: print *, "trying to assign ", int((prob(i) * numcells)+.499), & forrest@0: " cells for category", i forrest@0: jsum = jsum + int((prob(i) * numcells)+.499) forrest@0: end if forrest@0: 33 end do ! over all map categories forrest@0: ! forrest@0: print*, 'trying to assign a total of', jsum,' cells of', numcells, & forrest@0: ' total cells' forrest@0: !c forrest@0: !c set paint threshold for this map category forrest@0: ! hslice(i) = P(numcells - (int((prob(i) * numcells)+.499)),i) forrest@0: ! print *, "trying to assign ", int((prob(i) * numcells)+.499), & forrest@0: ! " cells for category", i forrest@0: ! end do ! over all map categories forrest@0: ! forrest@0: ! dump hslice array for checking forrest@0: print *, 'horizontal slices at these elevations:' forrest@0: do i = 0,ihabs forrest@0: print *, i, hslice(i) forrest@0: end do forrest@0: ! forrest@0: ! write XPM file of fractal probability terrain aspect forrest@0: if (probout) then forrest@0: do iz = 0, ihabs forrest@0: if (.not. (constraint(iz) .or. prob(iz) .eq. 0.0 .or. & forrest@0: iz .eq. igreatest)) then forrest@0: Call XPMRelief(P,iz,hslice(iz)) forrest@0: end if forrest@0: 335 end do forrest@0: end if forrest@0: ! forrest@0: ! forrest@0: numpainted = 0 forrest@0: numtied = 0 forrest@0: numblank = 0 forrest@0: ! forrest@0: print *,"Beginning tie resolution process" forrest@0: do l = 1, numcells ! start at top of P matrix forrest@0: idummy = 0 forrest@0: do j = 0, ihabs ! over all map categories forrest@0: ! probably don't need to do this if prob(j) .eq. 0.0 forrest@0: if (.not. (prob(j) .eq. 0.0 .or. j .eq. igreatest)) then forrest@0: ! forrest@0: ! if this cell is painted this category forrest@0: if ( P(l,j) .gt. hslice(j) ) then forrest@0: ! if ( P(l,j) .ge. hslice(j) ) then forrest@0: ! forrest@0: ! count the multiple assignments forrest@0: P(l,ihabs+3) = P(l,ihabs+3) + 1 forrest@0: ! forrest@0: ! find the winner category for this cell forrest@0: ! winner category has the highest probability surface forrest@0: if ( P(l,j) .gt. idummy ) then forrest@0: ihighest = j forrest@0: idummy = P(l,ihighest) forrest@0: end if forrest@0: ! forrest@0: end if ! this cell painted this category forrest@0: ! forrest@0: end if forrest@0: end do ! over all map categories forrest@0: ! print *, "highest for cell ", l," was category", ihighest forrest@0: ! forrest@0: ! keep count of each type of cell for displacement into forrest@0: ! sorted array later forrest@0: ! if singly-painted or currently tied, paint outmap with ihighest forrest@0: if ( P(l,ihabs+3) .eq. 1) then forrest@0: ! paint this cell with the winning category forrest@0: ! print *, "painting uncontested cell", P(l,ihabs+1), & forrest@0: ! P(l,ihabs+2)," with cat", ihighest forrest@0: ! forrest@0: numpainted = numpainted + 1 forrest@0: II( P(l,ihabs+1),P(l,ihabs+2) ) = ihighest forrest@0: else if ( P(l,ihabs+3) .eq. 0) then forrest@0: numblank = numblank + 1 forrest@0: else forrest@0: numtied = numtied + 1 forrest@0: ! print *, "painting tied cell", P(l,ihabs+1),P(l,ihabs+2), & forrest@0: ! " with winning cat", ihighest forrest@0: II( P(l,ihabs+1),P(l,ihabs+2) ) = ihighest forrest@0: end if forrest@0: ! forrest@0: end do forrest@0: ! forrest@0: ! at this point, II array has all winners and uncontested forrest@0: ! cells painted forrest@0: ! forrest@0: ! sort P array into 3 sections by number of assignments forrest@0: ! ascending sort results in empty, singly-painted, and forrest@0: ! ties sections forrest@0: print *,"Preparing to sort ties for resolution" forrest@0: ! forrest@0: ! sort by ties column forrest@0: istart = 1 forrest@0: Call HeapSort(P,numcells,istart,ihabs+3) forrest@0: ! print *, "got back from HeapSort" forrest@0: ! forrest@0: ! calculate pointers to section breaks forrest@0: mtsplit = numblank forrest@0: itiedsplit = numblank + numpainted forrest@0: ! forrest@0: mtresolved = mtsplit forrest@0: ! forrest@0: print *, "numblank is ", numblank forrest@0: print *, "numtied is ", numtied forrest@0: print *, "itiedsplit is ", itiedsplit forrest@0: print *, "numpainted is ", numpainted forrest@0: print *, "mtresolved and mtsplit are ", mtresolved forrest@0: forrest@0: if (numtied .gt. numblank) print*, 'WARNING: More tied cells', & forrest@0: ' than leftover blank cells to use for tie resolution!' forrest@0: ! forrest@0: ! sum multiwayness of ties and compare to number of empty cells forrest@0: do m = itiedsplit+1, numcells forrest@0: isum = isum + ( P(m,ihabs+3) - 1 ) forrest@0: end do forrest@0: print *, "sum of multiwayness of ties is ", isum forrest@0: if (isum .gt. numblank) print *, 'Warning: cells are tied more', & forrest@0: ' ways than there are blank cells to assign!' forrest@0: ! forrest@0: ! call dump_map(P) forrest@0: ! forrest@0: ! sort ties into priority order for resolution forrest@0: ! by calculating and sorting on summed elevation forrest@0: ! for all tied categories forrest@0: ! forrest@0: do m = itiedsplit+1, numcells ! over all tied forrest@0: kount = 0 forrest@0: idummy = 0 forrest@0: do k = 0, ihabs ! over all categories forrest@0: ! forrest@0: ! probably not necessary if prob(k) .eq. 0 forrest@0: if (.not. (prob(k) .eq. 0.0 .or. k .eq. igreatest)) then forrest@0: ! if this cell is painted this category forrest@0: if ( P(m,k) .gt. hslice(k) ) then forrest@0: ! if ( P(m,k) .ge. hslice(k) ) then forrest@0: kount = kount + 1 ! number of ties for this cell forrest@0: ties(kount) = k ! keep categories of ties forrest@0: ! forrest@0: end if forrest@0: end if forrest@0: 933 end do forrest@0: ! forrest@0: zmean = 0.0 forrest@0: do k = 1, kount forrest@0: zmean = zmean + P(m,ties(k)) forrest@0: end do forrest@0: ! P(m,ihabs+3) = zmean/kount forrest@0: ! try weighting priority by simple sum, not average forrest@0: ! therefore, max sum = (ihabs-1) * 32767 forrest@0: ! rescale as proportion of 32767 for integer*4 array cell forrest@0: P(m,ihabs+3) = zmean/(ihabs-1) forrest@0: ! forrest@0: end do forrest@0: ! forrest@0: ! sort ties by mean priority forrest@0: Call HeapSort(P,numtied,itiedsplit+1,ihabs+3) forrest@0: ! forrest@0: ! save min and max tie priorities for gray levels in XPMTies forrest@0: maxgray = P(numcells,ihabs+3) forrest@0: print*, "maxgray is ", maxgray forrest@0: ! forrest@0: ! generate xpm picture of tied cells forrest@0: ! colored by priority gray shades forrest@0: ! Call XPMTies(P) forrest@0: ! forrest@0: ! forrest@0: ! start at bottom of list (highest priority ties) forrest@0: do m = numcells, itiedsplit+1, -1 ! over all tied forrest@0: kount = 0 forrest@0: idummy = 0 forrest@0: do k = 0, ihabs ! over all categories forrest@0: ! forrest@0: ! probably not necessary if prob(k) .eq. 0 forrest@0: if (.not. (prob(k) .eq. 0.0 .or. k .eq. igreatest)) then forrest@0: ! if this cell is painted this category forrest@0: if ( P(m,k) .gt. hslice(k) ) then forrest@0: ! if ( P(m,k) .ge. hslice(k) ) then forrest@0: kount = kount + 1 ! number of ties for this cell forrest@0: ties(kount) = k ! keep categories of ties forrest@0: ! save most likely category for this cell forrest@0: if ( P(m,k) .gt. idummy ) then forrest@0: ihighest = k forrest@0: idummy = P(m,ihighest) forrest@0: end if forrest@0: end if forrest@0: end if forrest@0: 656 end do forrest@0: ! forrest@0: ! forrest@0: ! resolve all ties - forrest@0: ! over all tied categories that are not the winner forrest@0: ! (which has already been painted) forrest@0: ! we need an empty cell - forrest@0: ! search the empty cell list by this cat to find next forrest@0: ! most likely available empty cell forrest@0: ! forrest@0: ! empty cell list will be searched once for each tie category forrest@0: ! except the winning category forrest@0: ! forrest@0: ! print *, "deciding among ", kount, " ties for cell ", m forrest@0: do k = 1, kount ! over all ties for this cell forrest@0: if ( ties(k) .ne. ihighest ) then ! don't re-do the winner forrest@0: ! print *, ties(k), " was a tied category for cell ", m forrest@0: ! forrest@0: ibiggest = 0 forrest@0: idummy = 0 forrest@0: do i = 1,mtresolved forrest@0: if (P(i,ties(k)) .gt. idummy) then forrest@0: ibiggest = i forrest@0: idummy = P(i,ties(k)) forrest@0: end if forrest@0: end do forrest@0: ! forrest@0: Call Interchange(P,ibiggest,mtresolved) forrest@0: ! forrest@0: ! istart = 1 forrest@0: ! Call HeapSort(P,mtresolved,istart,ties(k)) forrest@0: ! forrest@0: ! now mtresolved points to highest forrest@0: ! next available empty cell of category ties(k) forrest@0: ! paint it! forrest@0: ! print *, "painting empty cell", mtresolved," at coords ", & forrest@0: ! P(mtresolved,ihabs+1), P(mtresolved,ihabs+2)," with prob ", & forrest@0: ! P(mtresolved,ties(k))," with cat", ties(k), & forrest@0: ! " for tie at cell", m forrest@0: II( P(mtresolved,ihabs+1),P(mtresolved,ihabs+2) ) = ties(k) forrest@0: ! burn one empty forrest@0: mtresolved = mtresolved - 1 forrest@0: end if forrest@0: ! forrest@0: 992 end do ! over all ties for this cell forrest@0: ! print *, "resolving ties at cell ", m forrest@0: ! forrest@0: end do ! over all tied cells forrest@0: ! forrest@0: !********************************************************************** forrest@0: ! forrest@0: ! forrest@0: ! Call dump_map(P) forrest@0: ! forrest@0: ! generate xpm picture of tied cells forrest@0: ! colored by priority gray shades forrest@0: ! and used blank cells colored in yellow forrest@0: Call XPMTies(P) forrest@0: ! forrest@0: ! forrest@0: ! set background category and forrest@0: ! perform summary of realized map forrest@0: do i = 0,N forrest@0: do j = 0,N forrest@0: if (II(i,j) .eq. 32767) II(i,j) = igreatest forrest@0: kcnt(II(i,j)) = kcnt(II(i,j)) + 1 forrest@0: ! write(12,*) II(i,j) forrest@0: end do forrest@0: end do forrest@0: ! forrest@0: ! if requested, forrest@0: ! write GRASS ascii map forrest@0: if (grassout) then forrest@0: write(10,91), N forrest@0: 91 format ('north: ',i4) forrest@0: write(10,92) forrest@0: 92 format ('south: 0') forrest@0: write(10,93), N forrest@0: 93 format ('east: ',i4) forrest@0: write(10,94) forrest@0: 94 format ('west: 0') forrest@0: write(10,95), N+1 forrest@0: 95 format ('rows: ',i4) forrest@0: write(10,96), N+1 forrest@0: 96 format ('cols: ',i4) forrest@0: write(fmt,'("(",I0,"(i5,1x))")') j forrest@0: do i = 0,N forrest@0: ! write(10,47) (II(i,j), j=0,N) forrest@0: write(10,fmt) (II(j,i), j=0,N) forrest@0: ! write(10,47) (II(j,i), j=0,N) forrest@0: !47 format ((i5,1x)) forrest@0: end do forrest@0: end if forrest@0: ! forrest@0: ! write final map to xpm file forrest@0: if (xpmout) Call XPMMap(N,II) forrest@0: ! forrest@0: print 299 forrest@0: 299 format (/t10,'Actual Map Summary:'/t12,'Habitat',5x, & forrest@0: 'P',3x,'Real Cells',3x,'Real P') forrest@0: ! pfract = 0.0 forrest@0: do i = 0, ihabs forrest@0: pfract = real(kcnt(i),kind(pfract)) / real(numcells,kind(pfract)) forrest@0: print 300, i, prob(i), kcnt(i), pfract forrest@0: 300 format (t10, i6, 4x, f6.5, 2x, i6, 2x, f6.5) forrest@0: end do forrest@0: ! forrest@0: !********************************************************************** forrest@0: ! forrest@0: ! visualize the synthetic landscape forrest@0: if (vizout) then forrest@0: ! forrest@0: ! print final map to screen forrest@0: do i=0,N forrest@0: do j=0,N forrest@0: idispvec(j) = II(i,j) forrest@0: end do forrest@0: ! print 199, (idispvec(j), j=0,N) forrest@0: !199 format ((i1, 1x)) forrest@0: end do forrest@0: ! forrest@0: ! print current state of P array forrest@0: ! Call print(P) forrest@0: ! forrest@0: ! visualize II array forrest@0: mapcolor = 'bobs.cmap\0' forrest@0: maptitle = 'Fractal Landscape Realizer \0' forrest@0: idtype = 4 forrest@0: idumm = 0 forrest@0: ipixid = -1 forrest@0: ! forrest@0: ! build display vector forrest@0: ! make small maps more visible by repeating cells and lines forrest@0: ! irepeat = maxN / N forrest@0: ! print *, irepeat forrest@0: k = 0 forrest@0: do i = 0,N forrest@0: ! kstart = k forrest@0: do j=0,N forrest@0: ! do kk = 1, irepeat forrest@0: ! idispvec(k)=II(i,j) forrest@0: idispvec(k)=II(j,i) forrest@0: k = k+1 forrest@0: ! end do ! repeating cells in this row forrest@0: end do ! all cells in this row forrest@0: ! kend = k-1 forrest@0: !c repeat lines forrest@0: ! if (irepeat .ge. 2) then forrest@0: ! do jj = 1,irepeat-1 forrest@0: ! do kk = kstart, kend forrest@0: ! idispvec(k) = idispvec(kk) forrest@0: ! k = k+1 forrest@0: ! end do ! dup as many times as necessary forrest@0: ! end do ! over all repeated lines forrest@0: ! end if forrest@0: ! forrest@0: end do forrest@0: ! forrest@0: !c adjust rows and columns by multiplier forrest@0: ! isize = (N+1) * irepeat forrest@0: ! jsize = (N+1) * irepeat forrest@0: ! print *, isize, jsize forrest@0: ! ierr=DrawPixmap(ipixid, isize, jsize, idtype, idispvec, forrest@0: ! ierr=DrawPixmap(ipixid,N+1,N+1,idtype,idispvec,idumm,mapcolor, & forrest@0: ! maptitle) forrest@0: ! print *,'ierr is', ierr forrest@0: ! print *, 'done visualizing' forrest@0: ! call sleep(500) forrest@0: ! forrest@0: end if ! if visualizing landscape forrest@0: ! forrest@0: stop forrest@0: ! end program MidPtFM2d forrest@0: end program realizer