1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/realizer_pro.f90 Wed Sep 26 17:50:53 2007 -0400
1.3 @@ -0,0 +1,942 @@
1.4 +program realizer
1.5 +! program MidPtFM2d
1.6 +!
1.7 +! fractal generator adapted from algorithm MidPointFM2D in the
1.8 +! Science of Fractal Images, M.F. Barnsley, R. L. Devaney
1.9 +! B. B. Mandelbrot, H.-O. Peitgen, D. Saupe, and R. F. Voss
1.10 +! Springer-Verlag
1.11 +!
1.12 +! this program modified to combine multiple fractal topographies
1.13 +! then normalize them relative to each other,
1.14 +! and perform multiple horizontal slices
1.15 +! to obtain requested p values for each habitat
1.16 +! in a new synthetic realization of the landscape
1.17 +!
1.18 +! William W. Hargrove
1.19 +! Paul M Schwartz
1.20 +! 5/7/96
1.21 +! (423) 241-2748
1.22 +!
1.23 + use realizer_mod
1.24 + implicit none
1.25 +
1.26 + integer(i4) :: P(maxN*maxN,0:maxsp+3)
1.27 + real(r4) :: X(0:maxN,0:maxN)
1.28 + real(r4) :: prob(0:maxsp)
1.29 + real(r4) :: ssum
1.30 + real(r4) :: ztemp
1.31 + integer(i4) :: igreatest
1.32 + integer(i4) :: hslice(0:maxsp)
1.33 + real(r4) :: elev(maxN*maxN)
1.34 + real(r4) :: H(maxN,maxsp)
1.35 + integer(i4) :: II(0:maxN,0:maxN)
1.36 + integer(i4) :: kcnt(0:maxsp)
1.37 + integer(i4) :: ties(maxsp)
1.38 + integer(i4) :: ifcnt
1.39 + integer(i4) :: iend
1.40 + integer(i4) :: istart
1.41 + integer(i4) :: icnt
1.42 + real(r4) :: sigma
1.43 + real(r4) :: xmin, xmax, xmean, xss
1.44 + real(r4) :: range1
1.45 + real(r4) :: zmean
1.46 + real(r4) :: pfract
1.47 + real(r4) :: xlow
1.48 + integer(i4) :: maxlevel
1.49 + integer(i4) :: numpainted
1.50 + integer(i4) :: numtied
1.51 + integer(i4) :: numblank
1.52 + integer(i4) :: ihighest
1.53 + integer(i4) :: ibiggest
1.54 + integer(i4) :: idummy
1.55 + integer(i4) :: idtype
1.56 + integer(i4) :: idumm
1.57 + integer(i4) :: isum, jsum
1.58 + integer(i4) :: ipixid
1.59 + integer(i4) :: kount
1.60 + integer(i4) :: i, j, k, l, m ! loop indices
1.61 + integer(i4) :: iv, iz
1.62 + integer(i4) :: N
1.63 + integer :: useed ! user seed value
1.64 + integer :: rsize ! length of array of seeds
1.65 + integer :: ierr ! error
1.66 + integer, allocatable :: iseed(:) ! array of random number seeds
1.67 + integer(i4) :: idispvec(0:maxN*maxN)
1.68 + character (len=1) :: ans
1.69 + character (len=28) :: maptitle
1.70 + character (len=60) :: mapfile
1.71 + character (len=10) :: mapcolor
1.72 + character (len=40) :: fmt ! For dynamic format strings in gfortran
1.73 + logical :: addition, wrap, gradient, slice, invert
1.74 + logical :: grassout, xpmout, probout, vizout
1.75 + logical :: constraint(0:maxsp)
1.76 +!
1.77 +! write execution script for next run
1.78 + open (9, file='input.scr',status='unknown')
1.79 +!
1.80 +!
1.81 +! enter random number seed
1.82 + print 196
1.83 +196 format (t10,'Do you want to enter a random number seed <Y>', &
1.84 + /t12,' or use the system clock? [clock]')
1.85 + read (5,97) ans
1.86 +97 format (a1)
1.87 + if (ans .eq. 'Y' .or. ans .eq. 'y') then
1.88 + write (9,197) ans
1.89 +197 format (a1,t20,"(enter seed or clock?)")
1.90 + print 100
1.91 +100 format (t10,'Enter random number seed ')
1.92 + read *, useed
1.93 + write (9,101) useed
1.94 +101 format (i12,t20,"(random number seed)")
1.95 + else
1.96 +21 useed = int(time8(),kind(useed))
1.97 + if (useed .eq. 0) then
1.98 + print *, "useed was zero."
1.99 + go to 21
1.100 + end if
1.101 + print 298
1.102 +298 format (t15,'Random number seed from system clock.'/)
1.103 + write (9,29) ans, useed
1.104 +29 format (a1,t20,"(random seed from clock = ",i12,")")
1.105 + end if
1.106 + ! Set the random number seed
1.107 + call random_seed(size=rsize)
1.108 + allocate(iseed(rsize),stat=ierr)
1.109 + if (ierr /= 0) then
1.110 + print *,'Error allocating iseed(rsize)'
1.111 + stop
1.112 + end if
1.113 + iseed(1:rsize) = useed
1.114 + call random_seed(put=iseed)
1.115 +!
1.116 + print 1017, useed
1.117 +1017 format (t15,'Random number seed = ',i12/)
1.118 +!
1.119 +!
1.120 + print 512
1.121 +512 format (t10,'Do you want to visualize', &
1.122 + /t12, ' the synthetic landscape? <Y,N>')
1.123 + read (5,522) ans
1.124 +522 format (a1)
1.125 + write (9,57) ans
1.126 +57 format (a1,t20,"(visualize the landscape?)")
1.127 + if (ans .eq. 'y' .or. ans .eq. 'Y') then
1.128 + vizout = .true.
1.129 + else
1.130 + vizout = .false.
1.131 + end if
1.132 +!
1.133 +1 print 225
1.134 +225 format (t10,'How many levels? (map is 2**levels on a side)')
1.135 + read *, maxlevel
1.136 +!
1.137 + N = 2**maxlevel
1.138 + irowcols = N
1.139 + numcells = (N+1)*(N+1)
1.140 + if (N .gt. maxN) goto 1
1.141 + write (9,102) maxlevel
1.142 +102 format (i3,t20,"(Number of levels)")
1.143 +!
1.144 + print 114
1.145 +114 format (t10,'Highest data category (excluding no-data)', &
1.146 + /t12,' in the realized synthetic landscape?')
1.147 + read *, ihabs
1.148 + write (9,199) ihabs
1.149 +199 format (i3,t20,"(highest data category)")
1.150 +!
1.151 + print 232
1.152 +232 format (t10,'Will you supply external probability maps', &
1.153 + /t12,' for any of these categories?')
1.154 + read (5,111) ans
1.155 +111 format (a1)
1.156 + write (9,161) ans
1.157 +161 format (a1,t20,"(supply external prob maps?)")
1.158 +!
1.159 + do j = 0, ihabs
1.160 + constraint(j) = .false.
1.161 + end do
1.162 +!
1.163 + isupply = 0
1.164 + if (ans .eq. 'Y' .or. ans .eq. 'y') then
1.165 +!
1.166 + print 134
1.167 +134 format (t10,'For how many categories in the synthetic', &
1.168 + ' landscape',/t12,' will you supply external probability', &
1.169 + ' maps?')
1.170 + read *, isupply
1.171 + write (9,189) isupply
1.172 +189 format (i3,t20,"(number of external prob maps)")
1.173 +!
1.174 + do i = 1, isupply
1.175 + print 216
1.176 +216 format (t10,'Constrain which map category?')
1.177 + read *, k
1.178 + write (9,799) k
1.179 +799 format (i3,t20,"(constrain which category?)")
1.180 +! set constraint flag for this category
1.181 + constraint(k) = .true.
1.182 +!
1.183 + print 220
1.184 +220 format (t10,'Enter name of input map file (integer*4, ', &
1.185 + /t12,' delimited by single spaces): ')
1.186 + read (5,230) mapfile
1.187 +230 format (a60)
1.188 + print 2030, mapfile
1.189 +2030 format (t15,'Input map name: ',a60/)
1.190 + write (9,231) mapfile
1.191 +231 format (a60)
1.192 +!
1.193 + print 144
1.194 +144 format (t10, 'Invert probabilities (i.e., DEM)? <Y,N>')
1.195 + read (5,321) ans
1.196 +321 format (a1)
1.197 + write (9,327) ans
1.198 +327 format (a1,t20,"(invert probabilities?)")
1.199 + if (ans .eq. 'y' .or. ans .eq. 'Y') then
1.200 + invert = .true.
1.201 + else
1.202 + invert = .false.
1.203 + end if
1.204 +!
1.205 +! load constraint maps in proper column of P array
1.206 +! invert probabilities if necessary, i.e., a DEM
1.207 + Call InputMap(mapfile,P,N,k,invert)
1.208 + end do
1.209 +!
1.210 + end if
1.211 +!
1.212 +!
1.213 + do iz = 0, ihabs
1.214 + if (.not. constraint(iz)) then
1.215 + do i = 1, maxlevel
1.216 +118 print 159, i, iz
1.217 +159 format (t10,'Input value of H for level ',i3,', category ',i3)
1.218 + read *, H(i,iz)
1.219 + print *, H(i,iz)
1.220 + write (9,155) H(i,iz), i, iz
1.221 +155 format (f4.2,t20,"(value of H for level ",i2,&
1.222 + ", category", i3,")")
1.223 + if (H(i,iz) .lt. 0.0 .or. H(i,iz) .gt. 1.0) goto 118
1.224 + end do
1.225 + end if
1.226 +664 end do
1.227 +!
1.228 +!
1.229 + ssum = 0.0
1.230 + do i = 1, ihabs
1.231 +115 print 116, i
1.232 +116 format (t10,'Probability for habitat ',i3)
1.233 + read *, prob(i)
1.234 + if (prob(i) .lt. 0.0) goto 115
1.235 + if (prob(i) .gt. 1.0) goto 115
1.236 + ssum = ssum + prob(i)
1.237 + if (ssum .gt. 1.0000) then
1.238 + print *, 'Total probability exceeds 1.0 - Re-enter'
1.239 + ssum = ssum - prob(i)
1.240 + goto 115
1.241 + end if
1.242 + write (9,151) prob(i), i
1.243 +151 format (f6.5,t20,"(probability of habitat type", i2,")")
1.244 + end do
1.245 +!
1.246 +! find greatest non-constraint p to make background
1.247 + ztemp = 0._r4
1.248 + igreatest = 0
1.249 + do i = 0, ihabs
1.250 + if (.not. constraint(i)) then
1.251 + if (prob(i) .gt. ztemp) then
1.252 + ztemp = prob(i)
1.253 + igreatest = i
1.254 + !print *, 'at i of ',i ,'new greatest is ', igreatest
1.255 + end if
1.256 + end if
1.257 +662 end do
1.258 +!
1.259 +! zero category by difference
1.260 + prob(0) = 1.0 - ssum
1.261 + if (prob(0) .lt. 0.00001) prob(0) = 0.0
1.262 + print *, 'prob(0) is ', prob(0)
1.263 +!
1.264 +! print *, 'igreatest is ', igreatest
1.265 + if (prob(0) .gt. prob(igreatest)) igreatest = 0
1.266 + print *, 'Cat ', igreatest,' will be calculated by absence'
1.267 +!
1.268 +! The Fractal Realizer assumes that the user-requested category
1.269 +! with the greatest coverage on the landscape (the highest requested p)
1.270 +! is the background or matrix category. In order to give the other
1.271 +! categories the greatest freedom for fractal tuning and to minimize
1.272 +! contention ties, this majority background color is computed by absence,
1.273 +! and all other categories are reclassed. Furthermore, categories which
1.274 +! are requested but for which the requested p = 0.0 are not actually
1.275 +! computed.
1.276 +!
1.277 +! calculate reclassed categories
1.278 +! ihabs = 0
1.279 +! do i = 0, ihabsmax
1.280 +! if (i .eq. igreatest) then
1.281 +! ireclass(0) = i
1.282 +! print *, 'Cat ', i,' reassigned to background'
1.283 +! goto 215
1.284 +! end if
1.285 +! if (prob(i) .eq. 0.0) then
1.286 +! print *, 'Cat ', i,' requested but not used - skipping'
1.287 +! goto 215
1.288 +! end if
1.289 +! ihabs = ihabs + 1
1.290 +! ireclass(ihabs) = i
1.291 +! print *, 'cat ', i,' will now be cat ', ihabs
1.292 +!
1.293 +!215 end do
1.294 +!
1.295 +!
1.296 +! do i = 0, ihabs
1.297 +! print 2002, ireclass(i), i, prob(ireclass(i))
1.298 +!2002 format ('Original category ', i2, ' re-mapped to ', &
1.299 +! i2, ' with requested probability of ', f6.4)
1.300 +! end do
1.301 +!
1.302 + print 212
1.303 +212 format (t10,'Generate GRASS r.in.ascii file? <Y,N>')
1.304 + read (5,222) ans
1.305 +222 format (a1)
1.306 + write (9,27) ans
1.307 +27 format (a1,t20,"(generate GRASS input file?)")
1.308 + if (ans .eq. 'y' .or. ans .eq. 'Y') then
1.309 + grassout = .true.
1.310 + else
1.311 + grassout = .false.
1.312 + end if
1.313 +!
1.314 + print 122
1.315 +122 format (t10,'Generate xpm file of synthetic landscape? <Y,N>')
1.316 + read (5,123) ans
1.317 +123 format (a1)
1.318 + write (9,124) ans
1.319 +124 format (a1,t20,"(xpm file of realized landscape?)")
1.320 + if (ans .eq. 'y' .or. ans .eq. 'Y') then
1.321 + xpmout = .true.
1.322 + else
1.323 + xpmout = .false.
1.324 + end if
1.325 +!
1.326 + print 322
1.327 +322 format (t10,'Generate xpm file of fractal probability' &
1.328 + ' landscapes? <Y,N>')
1.329 + read (5,323) ans
1.330 +323 format (a1)
1.331 + write (9,324) ans
1.332 +324 format (a1,t20,"(xpm file of probability landscapes?)")
1.333 + if (ans .eq. 'y' .or. ans .eq. 'Y') then
1.334 + probout = .true.
1.335 + else
1.336 + probout = .false.
1.337 + end if
1.338 +!
1.339 +! *****************************************************************
1.340 +! fill P array with a separate fractal probability map
1.341 +! for each desired habitat in the landscape
1.342 + do iz = 0, ihabs
1.343 + if (constraint(iz)) goto 22
1.344 +!
1.345 +! don't bother if igreatest or no cells requested of this category
1.346 + if (prob(iz) .eq. 0.0 .or. iz .eq. igreatest) then
1.347 + do i = 1, numcells
1.348 + P(i,iz) = 0
1.349 + end do
1.350 + goto 22
1.351 + end if
1.352 +!
1.353 +! print 102
1.354 +!102 format (t10,'Input the value for sigma')
1.355 +! read *, sigma
1.356 + sigma = 1.
1.357 + print 103
1.358 +103 format (t10,'Input the threshold for detection (xlow) ')
1.359 + read (5,119) xlow
1.360 +119 format (g12.6)
1.361 + write (9,125) xlow
1.362 +125 format (f12.6,t20,"(threshold for detection (xlow))")
1.363 + if (xlow .eq. 0) xlow = -1e30
1.364 + print *, xlow
1.365 + print 104
1.366 +104 format (t10,'Wrap? <Y,N>')
1.367 + read (5,145) ans
1.368 +145 format (a1)
1.369 + write (9,198) ans
1.370 +198 format (a1,t20,"(wrap?)")
1.371 + if (ans .eq. 'y' .or. ans .eq. 'Y') then
1.372 + wrap = .true.
1.373 + else
1.374 + wrap = .false.
1.375 + end if
1.376 +!
1.377 +! no gradient possible with wrap
1.378 + if (.not. wrap) then
1.379 + print 106
1.380 +106 format (t10,'Gradient? <Y,N>')
1.381 + read (5,107) ans
1.382 +107 format (a1)
1.383 + write (9,148) ans
1.384 +148 format (a1,t20,"(gradient?)")
1.385 + if (ans .eq. 'y' .or. ans .eq. 'Y') then
1.386 + gradient = .true.
1.387 +!
1.388 + print 1850
1.389 +1850 format (t10,'Gradient: Initial values for four map corners')
1.390 + print 1851
1.391 +1851 format (t10,'Northeast? ')
1.392 + read *, northeast
1.393 + print 1852, northeast
1.394 +1852 format (5x,g12.6)
1.395 + write (9,135) northeast
1.396 +135 format (f12.6,t20,"(Gradient: northeast)")
1.397 +!
1.398 + print 1853
1.399 +1853 format (t10,'Southeast? ')
1.400 + read *, southeast
1.401 + print 1852, southeast
1.402 + write (9,136) southeast
1.403 +136 format (f12.6,t20,"(Gradient: southeast)")
1.404 +!
1.405 + print 1854
1.406 +1854 format (t10,'Southwest? ')
1.407 + read *, southwest
1.408 + print 1852, southwest
1.409 + write (9,137) southwest
1.410 +137 format (f12.6,t20,"(Gradient: southwest)")
1.411 +!
1.412 + print 1855
1.413 +1855 format (t10,'Northwest? ')
1.414 + read *, northwest
1.415 + print 1852, northwest
1.416 + write (9,138) northwest
1.417 +138 format (f12.6,t20,"(Gradient: northwest)")
1.418 +!
1.419 + else
1.420 + gradient = .false.
1.421 + end if
1.422 +!
1.423 + end if !wrap = .false.
1.424 +!
1.425 +! print 111
1.426 +!111 format (t10,'Slice into discrete habitats? <Y,N>')
1.427 +! read (5,110) ans
1.428 +!110 format (a1)
1.429 +! if (ans .eq. 'y' .or. ans .eq. 'Y') then
1.430 +! slice = .true.
1.431 +! else
1.432 +! slice = .false.
1.433 +! end if
1.434 +!
1.435 + slice = .false.
1.436 +!
1.437 +! print 104
1.438 +! 104 format (t10,'Addition <Y,N>?')
1.439 +! read (5,105) ans
1.440 +! 105 format (a1)
1.441 +! if (ans .eq. 'y' .or. ans .eq. 'Y') then
1.442 +! addition = .true.
1.443 +! else
1.444 +! addition = .false.
1.445 +! end if
1.446 + addition = .false.
1.447 +!
1.448 +! call 2d multi-fractal generator
1.449 + print *,"Calling fractal generator for category ", iz
1.450 + call mpfm2d(X, maxlevel, sigma, H, addition, wrap, gradient, iz)
1.451 +!
1.452 + print *, "fractal generated"
1.453 + xmin = 1e30
1.454 + xmax = -1e30
1.455 + xmean = 0.0
1.456 + xss = 0.0
1.457 + icnt = 0
1.458 +!
1.459 +! obtain xmin and xmax for range normalization
1.460 +! calculate initial raw statistics
1.461 +! while you're at it
1.462 + do i = 0,N
1.463 + do j = 0,N
1.464 + if (X(i,j) .lt. xmin) xmin = X(i,j)
1.465 + if (X(i,j) .gt. xmax) xmax = X(i,j)
1.466 +! xmean = xmean + X(i,j)
1.467 +! xss = xss + X(i,j)*X(i,j)
1.468 +! icnt = icnt+1
1.469 + end do
1.470 + end do
1.471 +! xmean = xmean / icnt
1.472 +! print 200, xmean, xss, xmin, xmax
1.473 +!200 format (/t10,'Hfract Summary'/ &
1.474 +! t10,'Mean = ',g12.6 &
1.475 +! /t10,'Sum Sq. = ',g12.6 &
1.476 +! /t10,'Minimum = ',g12.6 &
1.477 +! /t10,'Maximum = ',g12.6)
1.478 +!
1.479 +! rescale and normalize fractal map between 0 and 32767, inclusive
1.480 +!
1.481 + range1 = xmax - xmin
1.482 + do i = 0,N
1.483 + do j = 0,N
1.484 + iv = (i*(N+1))+j+1
1.485 + P(iv,iz) = ((X(i,j) - xmin)/range1) * 32767
1.486 + P(iv,ihabs+1) = i
1.487 + P(iv,ihabs+2) = j
1.488 +! write(10,*) P(iv,iz)
1.489 + end do
1.490 + end do
1.491 + print*,'done normalizing now'
1.492 +!
1.493 +!c write GRASS ascii map
1.494 +! write(10,791), N
1.495 +!791 format ('north: ',i4)
1.496 +! write(10,792)
1.497 +!792 format ('south: 0')
1.498 +! write(10,793), N
1.499 +!793 format ('east: ',i4)
1.500 +! write(10,794)
1.501 +!794 format ('west: 0')
1.502 +! write(10,795), N+1
1.503 +!795 format ('rows: ',i4)
1.504 +! write(10,796), N+1
1.505 +!796 format ('cols: ',i4)
1.506 +! do i = 0,N
1.507 +! write(10,*) (P((i*(N+1))+j+1,iz), j=0,N)
1.508 +! end do
1.509 +!
1.510 +! end of iz loop over all landscape categories in synthetic map
1.511 +22 continue
1.512 + end do
1.513 +! **********************************************************************
1.514 +!
1.515 +! ihabs fractal maps now installed in P matrix:
1.516 +!
1.517 +! ihabs cols + 3
1.518 +! |---------------|
1.519 +! | |
1.520 +! | |
1.521 +! | empty cells |
1.522 +! mtresolved =>| |
1.523 +! | |
1.524 +! mtsplit =>|---------------|
1.525 +! | |
1.526 +! | |
1.527 +! | |
1.528 +! | |
1.529 +! | singly- |
1.530 +! | painted |
1.531 +! | |
1.532 +! | |
1.533 +! | |
1.534 +! | |
1.535 +! | |
1.536 +! itiedsplit =>|---------------|
1.537 +! | |
1.538 +! | |
1.539 +! | ties |
1.540 +! m =>| |
1.541 +! | |
1.542 +! numcells =>|_______________|
1.543 +!
1.544 +! ihabs+1 column carries x-coord of cell
1.545 +! ihabs+2 column carries y-coord of cell
1.546 +! ihabs+3 column carries number of ties for this cell
1.547 +!
1.548 +!
1.549 +! initialize II array with background value for re-use as output map
1.550 + II(:,:) = 32767
1.551 +!
1.552 +! set paint thresholds
1.553 + print *,"Calculating paint thresholds"
1.554 + do i = 0,ihabs ! over all map categories
1.555 +!
1.556 +!c don't bother if no cells requested of this cat
1.557 + if (prob(i) .ne. 0.0) then
1.558 +!
1.559 +! istart = 0
1.560 +! print *, "numcells is ", numcells
1.561 +! print *, "istart is ", istart
1.562 +! print *, "i is ", i
1.563 +! Call HeapSort(P,numcells,istart,i) ! sort P matrix by ith column
1.564 +!
1.565 +!
1.566 +! next category to be painted goes to single vector for sorting
1.567 + idispvec(0) = 0
1.568 + do l = 1, numcells
1.569 + idispvec(l) = P(l,i)
1.570 + end do
1.571 +!
1.572 +! call special vector version of Heapsort
1.573 + istart = 1
1.574 + Call HS2(idispvec,numcells+1,istart)
1.575 +!
1.576 +!c diagnostic print
1.577 +! do k=16636,16646
1.578 +! print *, i, k, idispvec(k)
1.579 +! end do
1.580 +!
1.581 +! set paint threshold for this map category
1.582 + hslice(i) = idispvec(numcells - (int((prob(i) * numcells)+.499)))
1.583 +! for absolute constraint maps
1.584 + if (hslice(i) .eq. 32767 .and. constraint(i)) hslice(i) = 32766
1.585 + print *, "trying to assign ", int((prob(i) * numcells)+.499), &
1.586 + " cells for category", i
1.587 + jsum = jsum + int((prob(i) * numcells)+.499)
1.588 + end if
1.589 +33 end do ! over all map categories
1.590 +!
1.591 + print*, 'trying to assign a total of', jsum,' cells of', numcells, &
1.592 + ' total cells'
1.593 +!c
1.594 +!c set paint threshold for this map category
1.595 +! hslice(i) = P(numcells - (int((prob(i) * numcells)+.499)),i)
1.596 +! print *, "trying to assign ", int((prob(i) * numcells)+.499), &
1.597 +! " cells for category", i
1.598 +! end do ! over all map categories
1.599 +!
1.600 +! dump hslice array for checking
1.601 + print *, 'horizontal slices at these elevations:'
1.602 + do i = 0,ihabs
1.603 + print *, i, hslice(i)
1.604 + end do
1.605 +!
1.606 +! write XPM file of fractal probability terrain aspect
1.607 + if (probout) then
1.608 + do iz = 0, ihabs
1.609 + if (.not. (constraint(iz) .or. prob(iz) .eq. 0.0 .or. &
1.610 + iz .eq. igreatest)) then
1.611 + Call XPMRelief(P,iz,hslice(iz))
1.612 + end if
1.613 +335 end do
1.614 + end if
1.615 +!
1.616 +!
1.617 + numpainted = 0
1.618 + numtied = 0
1.619 + numblank = 0
1.620 +!
1.621 + print *,"Beginning tie resolution process"
1.622 + do l = 1, numcells ! start at top of P matrix
1.623 + idummy = 0
1.624 + do j = 0, ihabs ! over all map categories
1.625 +! probably don't need to do this if prob(j) .eq. 0.0
1.626 + if (.not. (prob(j) .eq. 0.0 .or. j .eq. igreatest)) then
1.627 +!
1.628 +! if this cell is painted this category
1.629 + if ( P(l,j) .gt. hslice(j) ) then
1.630 +! if ( P(l,j) .ge. hslice(j) ) then
1.631 +!
1.632 +! count the multiple assignments
1.633 + P(l,ihabs+3) = P(l,ihabs+3) + 1
1.634 +!
1.635 +! find the winner category for this cell
1.636 +! winner category has the highest probability surface
1.637 + if ( P(l,j) .gt. idummy ) then
1.638 + ihighest = j
1.639 + idummy = P(l,ihighest)
1.640 + end if
1.641 +!
1.642 + end if ! this cell painted this category
1.643 +!
1.644 + end if
1.645 + end do ! over all map categories
1.646 +! print *, "highest for cell ", l," was category", ihighest
1.647 +!
1.648 +! keep count of each type of cell for displacement into
1.649 +! sorted array later
1.650 +! if singly-painted or currently tied, paint outmap with ihighest
1.651 + if ( P(l,ihabs+3) .eq. 1) then
1.652 +! paint this cell with the winning category
1.653 +! print *, "painting uncontested cell", P(l,ihabs+1), &
1.654 +! P(l,ihabs+2)," with cat", ihighest
1.655 +!
1.656 + numpainted = numpainted + 1
1.657 + II( P(l,ihabs+1),P(l,ihabs+2) ) = ihighest
1.658 + else if ( P(l,ihabs+3) .eq. 0) then
1.659 + numblank = numblank + 1
1.660 + else
1.661 + numtied = numtied + 1
1.662 +! print *, "painting tied cell", P(l,ihabs+1),P(l,ihabs+2), &
1.663 +! " with winning cat", ihighest
1.664 + II( P(l,ihabs+1),P(l,ihabs+2) ) = ihighest
1.665 + end if
1.666 +!
1.667 + end do
1.668 +!
1.669 +! at this point, II array has all winners and uncontested
1.670 +! cells painted
1.671 +!
1.672 +! sort P array into 3 sections by number of assignments
1.673 +! ascending sort results in empty, singly-painted, and
1.674 +! ties sections
1.675 + print *,"Preparing to sort ties for resolution"
1.676 +!
1.677 +! sort by ties column
1.678 + istart = 1
1.679 + Call HeapSort(P,numcells,istart,ihabs+3)
1.680 +! print *, "got back from HeapSort"
1.681 +!
1.682 +! calculate pointers to section breaks
1.683 + mtsplit = numblank
1.684 + itiedsplit = numblank + numpainted
1.685 +!
1.686 + mtresolved = mtsplit
1.687 +!
1.688 + print *, "numblank is ", numblank
1.689 + print *, "numtied is ", numtied
1.690 + print *, "itiedsplit is ", itiedsplit
1.691 + print *, "numpainted is ", numpainted
1.692 + print *, "mtresolved and mtsplit are ", mtresolved
1.693 +
1.694 + if (numtied .gt. numblank) print*, 'WARNING: More tied cells', &
1.695 + ' than leftover blank cells to use for tie resolution!'
1.696 +!
1.697 +! sum multiwayness of ties and compare to number of empty cells
1.698 + do m = itiedsplit+1, numcells
1.699 + isum = isum + ( P(m,ihabs+3) - 1 )
1.700 + end do
1.701 + print *, "sum of multiwayness of ties is ", isum
1.702 + if (isum .gt. numblank) print *, 'Warning: cells are tied more', &
1.703 + ' ways than there are blank cells to assign!'
1.704 +!
1.705 +! call dump_map(P)
1.706 +!
1.707 +! sort ties into priority order for resolution
1.708 +! by calculating and sorting on summed elevation
1.709 +! for all tied categories
1.710 +!
1.711 + do m = itiedsplit+1, numcells ! over all tied
1.712 + kount = 0
1.713 + idummy = 0
1.714 + do k = 0, ihabs ! over all categories
1.715 +!
1.716 +! probably not necessary if prob(k) .eq. 0
1.717 + if (.not. (prob(k) .eq. 0.0 .or. k .eq. igreatest)) then
1.718 +! if this cell is painted this category
1.719 + if ( P(m,k) .gt. hslice(k) ) then
1.720 +! if ( P(m,k) .ge. hslice(k) ) then
1.721 + kount = kount + 1 ! number of ties for this cell
1.722 + ties(kount) = k ! keep categories of ties
1.723 +!
1.724 + end if
1.725 + end if
1.726 +933 end do
1.727 +!
1.728 + zmean = 0.0
1.729 + do k = 1, kount
1.730 + zmean = zmean + P(m,ties(k))
1.731 + end do
1.732 +! P(m,ihabs+3) = zmean/kount
1.733 +! try weighting priority by simple sum, not average
1.734 +! therefore, max sum = (ihabs-1) * 32767
1.735 +! rescale as proportion of 32767 for integer*4 array cell
1.736 + P(m,ihabs+3) = zmean/(ihabs-1)
1.737 +!
1.738 + end do
1.739 +!
1.740 +! sort ties by mean priority
1.741 + Call HeapSort(P,numtied,itiedsplit+1,ihabs+3)
1.742 +!
1.743 +! save min and max tie priorities for gray levels in XPMTies
1.744 + maxgray = P(numcells,ihabs+3)
1.745 + print*, "maxgray is ", maxgray
1.746 +!
1.747 +! generate xpm picture of tied cells
1.748 +! colored by priority gray shades
1.749 +! Call XPMTies(P)
1.750 +!
1.751 +!
1.752 +! start at bottom of list (highest priority ties)
1.753 + do m = numcells, itiedsplit+1, -1 ! over all tied
1.754 + kount = 0
1.755 + idummy = 0
1.756 + do k = 0, ihabs ! over all categories
1.757 +!
1.758 +! probably not necessary if prob(k) .eq. 0
1.759 + if (.not. (prob(k) .eq. 0.0 .or. k .eq. igreatest)) then
1.760 +! if this cell is painted this category
1.761 + if ( P(m,k) .gt. hslice(k) ) then
1.762 +! if ( P(m,k) .ge. hslice(k) ) then
1.763 + kount = kount + 1 ! number of ties for this cell
1.764 + ties(kount) = k ! keep categories of ties
1.765 +! save most likely category for this cell
1.766 + if ( P(m,k) .gt. idummy ) then
1.767 + ihighest = k
1.768 + idummy = P(m,ihighest)
1.769 + end if
1.770 + end if
1.771 + end if
1.772 +656 end do
1.773 +!
1.774 +!
1.775 +! resolve all ties -
1.776 +! over all tied categories that are not the winner
1.777 +! (which has already been painted)
1.778 +! we need an empty cell -
1.779 +! search the empty cell list by this cat to find next
1.780 +! most likely available empty cell
1.781 +!
1.782 +! empty cell list will be searched once for each tie category
1.783 +! except the winning category
1.784 +!
1.785 +! print *, "deciding among ", kount, " ties for cell ", m
1.786 + do k = 1, kount ! over all ties for this cell
1.787 + if ( ties(k) .ne. ihighest ) then ! don't re-do the winner
1.788 +! print *, ties(k), " was a tied category for cell ", m
1.789 +!
1.790 + ibiggest = 0
1.791 + idummy = 0
1.792 + do i = 1,mtresolved
1.793 + if (P(i,ties(k)) .gt. idummy) then
1.794 + ibiggest = i
1.795 + idummy = P(i,ties(k))
1.796 + end if
1.797 + end do
1.798 +!
1.799 + Call Interchange(P,ibiggest,mtresolved)
1.800 +!
1.801 +! istart = 1
1.802 +! Call HeapSort(P,mtresolved,istart,ties(k))
1.803 +!
1.804 +! now mtresolved points to highest
1.805 +! next available empty cell of category ties(k)
1.806 +! paint it!
1.807 +! print *, "painting empty cell", mtresolved," at coords ", &
1.808 +! P(mtresolved,ihabs+1), P(mtresolved,ihabs+2)," with prob ", &
1.809 +! P(mtresolved,ties(k))," with cat", ties(k), &
1.810 +! " for tie at cell", m
1.811 + II( P(mtresolved,ihabs+1),P(mtresolved,ihabs+2) ) = ties(k)
1.812 +! burn one empty
1.813 + mtresolved = mtresolved - 1
1.814 + end if
1.815 +!
1.816 +992 end do ! over all ties for this cell
1.817 +! print *, "resolving ties at cell ", m
1.818 +!
1.819 + end do ! over all tied cells
1.820 +!
1.821 +!**********************************************************************
1.822 +!
1.823 +!
1.824 +! Call dump_map(P)
1.825 +!
1.826 +! generate xpm picture of tied cells
1.827 +! colored by priority gray shades
1.828 +! and used blank cells colored in yellow
1.829 + Call XPMTies(P)
1.830 +!
1.831 +!
1.832 +! set background category and
1.833 +! perform summary of realized map
1.834 + do i = 0,N
1.835 + do j = 0,N
1.836 + if (II(i,j) .eq. 32767) II(i,j) = igreatest
1.837 + kcnt(II(i,j)) = kcnt(II(i,j)) + 1
1.838 +! write(12,*) II(i,j)
1.839 + end do
1.840 + end do
1.841 +!
1.842 +! if requested,
1.843 +! write GRASS ascii map
1.844 + if (grassout) then
1.845 + write(10,91), N
1.846 +91 format ('north: ',i4)
1.847 + write(10,92)
1.848 +92 format ('south: 0')
1.849 + write(10,93), N
1.850 +93 format ('east: ',i4)
1.851 + write(10,94)
1.852 +94 format ('west: 0')
1.853 + write(10,95), N+1
1.854 +95 format ('rows: ',i4)
1.855 + write(10,96), N+1
1.856 +96 format ('cols: ',i4)
1.857 + write(fmt,'("(",I0,"(i5,1x))")') j
1.858 + do i = 0,N
1.859 +! write(10,47) (II(i,j), j=0,N)
1.860 + write(10,fmt) (II(j,i), j=0,N)
1.861 +! write(10,47) (II(j,i), j=0,N)
1.862 +!47 format (<j>(i5,1x))
1.863 + end do
1.864 + end if
1.865 +!
1.866 +! write final map to xpm file
1.867 + if (xpmout) Call XPMMap(N,II)
1.868 +!
1.869 + print 299
1.870 +299 format (/t10,'Actual Map Summary:'/t12,'Habitat',5x, &
1.871 + 'P',3x,'Real Cells',3x,'Real P')
1.872 +! pfract = 0.0
1.873 + do i = 0, ihabs
1.874 + pfract = real(kcnt(i),kind(pfract)) / real(numcells,kind(pfract))
1.875 + print 300, i, prob(i), kcnt(i), pfract
1.876 +300 format (t10, i6, 4x, f6.5, 2x, i6, 2x, f6.5)
1.877 + end do
1.878 +!
1.879 +!**********************************************************************
1.880 +!
1.881 +! visualize the synthetic landscape
1.882 + if (vizout) then
1.883 +!
1.884 +! print final map to screen
1.885 + do i=0,N
1.886 + do j=0,N
1.887 + idispvec(j) = II(i,j)
1.888 + end do
1.889 +! print 199, (idispvec(j), j=0,N)
1.890 +!199 format (<N+1>(i1, 1x))
1.891 + end do
1.892 +!
1.893 +! print current state of P array
1.894 +! Call print(P)
1.895 +!
1.896 +! visualize II array
1.897 + mapcolor = 'bobs.cmap\0'
1.898 + maptitle = 'Fractal Landscape Realizer \0'
1.899 + idtype = 4
1.900 + idumm = 0
1.901 + ipixid = -1
1.902 +!
1.903 +! build display vector
1.904 +! make small maps more visible by repeating cells and lines
1.905 +! irepeat = maxN / N
1.906 +! print *, irepeat
1.907 + k = 0
1.908 + do i = 0,N
1.909 +! kstart = k
1.910 + do j=0,N
1.911 +! do kk = 1, irepeat
1.912 +! idispvec(k)=II(i,j)
1.913 + idispvec(k)=II(j,i)
1.914 + k = k+1
1.915 +! end do ! repeating cells in this row
1.916 + end do ! all cells in this row
1.917 +! kend = k-1
1.918 +!c repeat lines
1.919 +! if (irepeat .ge. 2) then
1.920 +! do jj = 1,irepeat-1
1.921 +! do kk = kstart, kend
1.922 +! idispvec(k) = idispvec(kk)
1.923 +! k = k+1
1.924 +! end do ! dup as many times as necessary
1.925 +! end do ! over all repeated lines
1.926 +! end if
1.927 +!
1.928 + end do
1.929 +!
1.930 +!c adjust rows and columns by multiplier
1.931 +! isize = (N+1) * irepeat
1.932 +! jsize = (N+1) * irepeat
1.933 +! print *, isize, jsize
1.934 +! ierr=DrawPixmap(ipixid, isize, jsize, idtype, idispvec,
1.935 +! ierr=DrawPixmap(ipixid,N+1,N+1,idtype,idispvec,idumm,mapcolor, &
1.936 +! maptitle)
1.937 +! print *,'ierr is', ierr
1.938 +! print *, 'done visualizing'
1.939 +! call sleep(500)
1.940 +!
1.941 + end if ! if visualizing landscape
1.942 +!
1.943 + stop
1.944 +! end program MidPtFM2d
1.945 +end program realizer