|
1 program realizer |
|
2 ! program MidPtFM2d |
|
3 ! |
|
4 ! fractal generator adapted from algorithm MidPointFM2D in the |
|
5 ! Science of Fractal Images, M.F. Barnsley, R. L. Devaney |
|
6 ! B. B. Mandelbrot, H.-O. Peitgen, D. Saupe, and R. F. Voss |
|
7 ! Springer-Verlag |
|
8 ! |
|
9 ! this program modified to combine multiple fractal topographies |
|
10 ! then normalize them relative to each other, |
|
11 ! and perform multiple horizontal slices |
|
12 ! to obtain requested p values for each habitat |
|
13 ! in a new synthetic realization of the landscape |
|
14 ! |
|
15 ! William W. Hargrove |
|
16 ! Paul M Schwartz |
|
17 ! 5/7/96 |
|
18 ! (423) 241-2748 |
|
19 ! |
|
20 use realizer_mod |
|
21 implicit none |
|
22 |
|
23 integer(i4) :: P(maxN*maxN,0:maxsp+3) |
|
24 real(r4) :: X(0:maxN,0:maxN) |
|
25 real(r4) :: prob(0:maxsp) |
|
26 real(r4) :: ssum |
|
27 real(r4) :: ztemp |
|
28 integer(i4) :: igreatest |
|
29 integer(i4) :: hslice(0:maxsp) |
|
30 real(r4) :: elev(maxN*maxN) |
|
31 real(r4) :: H(maxN,maxsp) |
|
32 integer(i4) :: II(0:maxN,0:maxN) |
|
33 integer(i4) :: kcnt(0:maxsp) |
|
34 integer(i4) :: ties(maxsp) |
|
35 integer(i4) :: ifcnt |
|
36 integer(i4) :: iend |
|
37 integer(i4) :: istart |
|
38 integer(i4) :: icnt |
|
39 real(r4) :: sigma |
|
40 real(r4) :: xmin, xmax, xmean, xss |
|
41 real(r4) :: range1 |
|
42 real(r4) :: zmean |
|
43 real(r4) :: pfract |
|
44 real(r4) :: xlow |
|
45 integer(i4) :: maxlevel |
|
46 integer(i4) :: numpainted |
|
47 integer(i4) :: numtied |
|
48 integer(i4) :: numblank |
|
49 integer(i4) :: ihighest |
|
50 integer(i4) :: ibiggest |
|
51 integer(i4) :: idummy |
|
52 integer(i4) :: idtype |
|
53 integer(i4) :: idumm |
|
54 integer(i4) :: isum, jsum |
|
55 integer(i4) :: ipixid |
|
56 integer(i4) :: kount |
|
57 integer(i4) :: i, j, k, l, m ! loop indices |
|
58 integer(i4) :: iv, iz |
|
59 integer(i4) :: N |
|
60 integer :: useed ! user seed value |
|
61 integer :: rsize ! length of array of seeds |
|
62 integer :: ierr ! error |
|
63 integer, allocatable :: iseed(:) ! array of random number seeds |
|
64 integer(i4) :: idispvec(0:maxN*maxN) |
|
65 character (len=1) :: ans |
|
66 character (len=28) :: maptitle |
|
67 character (len=60) :: mapfile |
|
68 character (len=10) :: mapcolor |
|
69 character (len=40) :: fmt ! For dynamic format strings in gfortran |
|
70 logical :: addition, wrap, gradient, slice, invert |
|
71 logical :: grassout, xpmout, probout, vizout |
|
72 logical :: constraint(0:maxsp) |
|
73 ! |
|
74 ! write execution script for next run |
|
75 open (9, file='input.scr',status='unknown') |
|
76 ! |
|
77 ! |
|
78 ! enter random number seed |
|
79 print 196 |
|
80 196 format (t10,'Do you want to enter a random number seed <Y>', & |
|
81 /t12,' or use the system clock? [clock]') |
|
82 read (5,97) ans |
|
83 97 format (a1) |
|
84 if (ans .eq. 'Y' .or. ans .eq. 'y') then |
|
85 write (9,197) ans |
|
86 197 format (a1,t20,"(enter seed or clock?)") |
|
87 print 100 |
|
88 100 format (t10,'Enter random number seed ') |
|
89 read *, useed |
|
90 write (9,101) useed |
|
91 101 format (i12,t20,"(random number seed)") |
|
92 else |
|
93 21 useed = int(time8(),kind(useed)) |
|
94 if (useed .eq. 0) then |
|
95 print *, "useed was zero." |
|
96 go to 21 |
|
97 end if |
|
98 print 298 |
|
99 298 format (t15,'Random number seed from system clock.'/) |
|
100 write (9,29) ans, useed |
|
101 29 format (a1,t20,"(random seed from clock = ",i12,")") |
|
102 end if |
|
103 ! Set the random number seed |
|
104 call random_seed(size=rsize) |
|
105 allocate(iseed(rsize),stat=ierr) |
|
106 if (ierr /= 0) then |
|
107 print *,'Error allocating iseed(rsize)' |
|
108 stop |
|
109 end if |
|
110 iseed(1:rsize) = useed |
|
111 call random_seed(put=iseed) |
|
112 ! |
|
113 print 1017, useed |
|
114 1017 format (t15,'Random number seed = ',i12/) |
|
115 ! |
|
116 ! |
|
117 print 512 |
|
118 512 format (t10,'Do you want to visualize', & |
|
119 /t12, ' the synthetic landscape? <Y,N>') |
|
120 read (5,522) ans |
|
121 522 format (a1) |
|
122 write (9,57) ans |
|
123 57 format (a1,t20,"(visualize the landscape?)") |
|
124 if (ans .eq. 'y' .or. ans .eq. 'Y') then |
|
125 vizout = .true. |
|
126 else |
|
127 vizout = .false. |
|
128 end if |
|
129 ! |
|
130 1 print 225 |
|
131 225 format (t10,'How many levels? (map is 2**levels on a side)') |
|
132 read *, maxlevel |
|
133 ! |
|
134 N = 2**maxlevel |
|
135 irowcols = N |
|
136 numcells = (N+1)*(N+1) |
|
137 if (N .gt. maxN) goto 1 |
|
138 write (9,102) maxlevel |
|
139 102 format (i3,t20,"(Number of levels)") |
|
140 ! |
|
141 print 114 |
|
142 114 format (t10,'Highest data category (excluding no-data)', & |
|
143 /t12,' in the realized synthetic landscape?') |
|
144 read *, ihabs |
|
145 write (9,199) ihabs |
|
146 199 format (i3,t20,"(highest data category)") |
|
147 ! |
|
148 print 232 |
|
149 232 format (t10,'Will you supply external probability maps', & |
|
150 /t12,' for any of these categories?') |
|
151 read (5,111) ans |
|
152 111 format (a1) |
|
153 write (9,161) ans |
|
154 161 format (a1,t20,"(supply external prob maps?)") |
|
155 ! |
|
156 do j = 0, ihabs |
|
157 constraint(j) = .false. |
|
158 end do |
|
159 ! |
|
160 isupply = 0 |
|
161 if (ans .eq. 'Y' .or. ans .eq. 'y') then |
|
162 ! |
|
163 print 134 |
|
164 134 format (t10,'For how many categories in the synthetic', & |
|
165 ' landscape',/t12,' will you supply external probability', & |
|
166 ' maps?') |
|
167 read *, isupply |
|
168 write (9,189) isupply |
|
169 189 format (i3,t20,"(number of external prob maps)") |
|
170 ! |
|
171 do i = 1, isupply |
|
172 print 216 |
|
173 216 format (t10,'Constrain which map category?') |
|
174 read *, k |
|
175 write (9,799) k |
|
176 799 format (i3,t20,"(constrain which category?)") |
|
177 ! set constraint flag for this category |
|
178 constraint(k) = .true. |
|
179 ! |
|
180 print 220 |
|
181 220 format (t10,'Enter name of input map file (integer*4, ', & |
|
182 /t12,' delimited by single spaces): ') |
|
183 read (5,230) mapfile |
|
184 230 format (a60) |
|
185 print 2030, mapfile |
|
186 2030 format (t15,'Input map name: ',a60/) |
|
187 write (9,231) mapfile |
|
188 231 format (a60) |
|
189 ! |
|
190 print 144 |
|
191 144 format (t10, 'Invert probabilities (i.e., DEM)? <Y,N>') |
|
192 read (5,321) ans |
|
193 321 format (a1) |
|
194 write (9,327) ans |
|
195 327 format (a1,t20,"(invert probabilities?)") |
|
196 if (ans .eq. 'y' .or. ans .eq. 'Y') then |
|
197 invert = .true. |
|
198 else |
|
199 invert = .false. |
|
200 end if |
|
201 ! |
|
202 ! load constraint maps in proper column of P array |
|
203 ! invert probabilities if necessary, i.e., a DEM |
|
204 Call InputMap(mapfile,P,N,k,invert) |
|
205 end do |
|
206 ! |
|
207 end if |
|
208 ! |
|
209 ! |
|
210 do iz = 0, ihabs |
|
211 if (.not. constraint(iz)) then |
|
212 do i = 1, maxlevel |
|
213 118 print 159, i, iz |
|
214 159 format (t10,'Input value of H for level ',i3,', category ',i3) |
|
215 read *, H(i,iz) |
|
216 print *, H(i,iz) |
|
217 write (9,155) H(i,iz), i, iz |
|
218 155 format (f4.2,t20,"(value of H for level ",i2,& |
|
219 ", category", i3,")") |
|
220 if (H(i,iz) .lt. 0.0 .or. H(i,iz) .gt. 1.0) goto 118 |
|
221 end do |
|
222 end if |
|
223 664 end do |
|
224 ! |
|
225 ! |
|
226 ssum = 0.0 |
|
227 do i = 1, ihabs |
|
228 115 print 116, i |
|
229 116 format (t10,'Probability for habitat ',i3) |
|
230 read *, prob(i) |
|
231 if (prob(i) .lt. 0.0) goto 115 |
|
232 if (prob(i) .gt. 1.0) goto 115 |
|
233 ssum = ssum + prob(i) |
|
234 if (ssum .gt. 1.0000) then |
|
235 print *, 'Total probability exceeds 1.0 - Re-enter' |
|
236 ssum = ssum - prob(i) |
|
237 goto 115 |
|
238 end if |
|
239 write (9,151) prob(i), i |
|
240 151 format (f6.5,t20,"(probability of habitat type", i2,")") |
|
241 end do |
|
242 ! |
|
243 ! find greatest non-constraint p to make background |
|
244 ztemp = 0._r4 |
|
245 igreatest = 0 |
|
246 do i = 0, ihabs |
|
247 if (.not. constraint(i)) then |
|
248 if (prob(i) .gt. ztemp) then |
|
249 ztemp = prob(i) |
|
250 igreatest = i |
|
251 !print *, 'at i of ',i ,'new greatest is ', igreatest |
|
252 end if |
|
253 end if |
|
254 662 end do |
|
255 ! |
|
256 ! zero category by difference |
|
257 prob(0) = 1.0 - ssum |
|
258 if (prob(0) .lt. 0.00001) prob(0) = 0.0 |
|
259 print *, 'prob(0) is ', prob(0) |
|
260 ! |
|
261 ! print *, 'igreatest is ', igreatest |
|
262 if (prob(0) .gt. prob(igreatest)) igreatest = 0 |
|
263 print *, 'Cat ', igreatest,' will be calculated by absence' |
|
264 ! |
|
265 ! The Fractal Realizer assumes that the user-requested category |
|
266 ! with the greatest coverage on the landscape (the highest requested p) |
|
267 ! is the background or matrix category. In order to give the other |
|
268 ! categories the greatest freedom for fractal tuning and to minimize |
|
269 ! contention ties, this majority background color is computed by absence, |
|
270 ! and all other categories are reclassed. Furthermore, categories which |
|
271 ! are requested but for which the requested p = 0.0 are not actually |
|
272 ! computed. |
|
273 ! |
|
274 ! calculate reclassed categories |
|
275 ! ihabs = 0 |
|
276 ! do i = 0, ihabsmax |
|
277 ! if (i .eq. igreatest) then |
|
278 ! ireclass(0) = i |
|
279 ! print *, 'Cat ', i,' reassigned to background' |
|
280 ! goto 215 |
|
281 ! end if |
|
282 ! if (prob(i) .eq. 0.0) then |
|
283 ! print *, 'Cat ', i,' requested but not used - skipping' |
|
284 ! goto 215 |
|
285 ! end if |
|
286 ! ihabs = ihabs + 1 |
|
287 ! ireclass(ihabs) = i |
|
288 ! print *, 'cat ', i,' will now be cat ', ihabs |
|
289 ! |
|
290 !215 end do |
|
291 ! |
|
292 ! |
|
293 ! do i = 0, ihabs |
|
294 ! print 2002, ireclass(i), i, prob(ireclass(i)) |
|
295 !2002 format ('Original category ', i2, ' re-mapped to ', & |
|
296 ! i2, ' with requested probability of ', f6.4) |
|
297 ! end do |
|
298 ! |
|
299 print 212 |
|
300 212 format (t10,'Generate GRASS r.in.ascii file? <Y,N>') |
|
301 read (5,222) ans |
|
302 222 format (a1) |
|
303 write (9,27) ans |
|
304 27 format (a1,t20,"(generate GRASS input file?)") |
|
305 if (ans .eq. 'y' .or. ans .eq. 'Y') then |
|
306 grassout = .true. |
|
307 else |
|
308 grassout = .false. |
|
309 end if |
|
310 ! |
|
311 print 122 |
|
312 122 format (t10,'Generate xpm file of synthetic landscape? <Y,N>') |
|
313 read (5,123) ans |
|
314 123 format (a1) |
|
315 write (9,124) ans |
|
316 124 format (a1,t20,"(xpm file of realized landscape?)") |
|
317 if (ans .eq. 'y' .or. ans .eq. 'Y') then |
|
318 xpmout = .true. |
|
319 else |
|
320 xpmout = .false. |
|
321 end if |
|
322 ! |
|
323 print 322 |
|
324 322 format (t10,'Generate xpm file of fractal probability' & |
|
325 ' landscapes? <Y,N>') |
|
326 read (5,323) ans |
|
327 323 format (a1) |
|
328 write (9,324) ans |
|
329 324 format (a1,t20,"(xpm file of probability landscapes?)") |
|
330 if (ans .eq. 'y' .or. ans .eq. 'Y') then |
|
331 probout = .true. |
|
332 else |
|
333 probout = .false. |
|
334 end if |
|
335 ! |
|
336 ! ***************************************************************** |
|
337 ! fill P array with a separate fractal probability map |
|
338 ! for each desired habitat in the landscape |
|
339 do iz = 0, ihabs |
|
340 if (constraint(iz)) goto 22 |
|
341 ! |
|
342 ! don't bother if igreatest or no cells requested of this category |
|
343 if (prob(iz) .eq. 0.0 .or. iz .eq. igreatest) then |
|
344 do i = 1, numcells |
|
345 P(i,iz) = 0 |
|
346 end do |
|
347 goto 22 |
|
348 end if |
|
349 ! |
|
350 ! print 102 |
|
351 !102 format (t10,'Input the value for sigma') |
|
352 ! read *, sigma |
|
353 sigma = 1. |
|
354 print 103 |
|
355 103 format (t10,'Input the threshold for detection (xlow) ') |
|
356 read (5,119) xlow |
|
357 119 format (g12.6) |
|
358 write (9,125) xlow |
|
359 125 format (f12.6,t20,"(threshold for detection (xlow))") |
|
360 if (xlow .eq. 0) xlow = -1e30 |
|
361 print *, xlow |
|
362 print 104 |
|
363 104 format (t10,'Wrap? <Y,N>') |
|
364 read (5,145) ans |
|
365 145 format (a1) |
|
366 write (9,198) ans |
|
367 198 format (a1,t20,"(wrap?)") |
|
368 if (ans .eq. 'y' .or. ans .eq. 'Y') then |
|
369 wrap = .true. |
|
370 else |
|
371 wrap = .false. |
|
372 end if |
|
373 ! |
|
374 ! no gradient possible with wrap |
|
375 if (.not. wrap) then |
|
376 print 106 |
|
377 106 format (t10,'Gradient? <Y,N>') |
|
378 read (5,107) ans |
|
379 107 format (a1) |
|
380 write (9,148) ans |
|
381 148 format (a1,t20,"(gradient?)") |
|
382 if (ans .eq. 'y' .or. ans .eq. 'Y') then |
|
383 gradient = .true. |
|
384 ! |
|
385 print 1850 |
|
386 1850 format (t10,'Gradient: Initial values for four map corners') |
|
387 print 1851 |
|
388 1851 format (t10,'Northeast? ') |
|
389 read *, northeast |
|
390 print 1852, northeast |
|
391 1852 format (5x,g12.6) |
|
392 write (9,135) northeast |
|
393 135 format (f12.6,t20,"(Gradient: northeast)") |
|
394 ! |
|
395 print 1853 |
|
396 1853 format (t10,'Southeast? ') |
|
397 read *, southeast |
|
398 print 1852, southeast |
|
399 write (9,136) southeast |
|
400 136 format (f12.6,t20,"(Gradient: southeast)") |
|
401 ! |
|
402 print 1854 |
|
403 1854 format (t10,'Southwest? ') |
|
404 read *, southwest |
|
405 print 1852, southwest |
|
406 write (9,137) southwest |
|
407 137 format (f12.6,t20,"(Gradient: southwest)") |
|
408 ! |
|
409 print 1855 |
|
410 1855 format (t10,'Northwest? ') |
|
411 read *, northwest |
|
412 print 1852, northwest |
|
413 write (9,138) northwest |
|
414 138 format (f12.6,t20,"(Gradient: northwest)") |
|
415 ! |
|
416 else |
|
417 gradient = .false. |
|
418 end if |
|
419 ! |
|
420 end if !wrap = .false. |
|
421 ! |
|
422 ! print 111 |
|
423 !111 format (t10,'Slice into discrete habitats? <Y,N>') |
|
424 ! read (5,110) ans |
|
425 !110 format (a1) |
|
426 ! if (ans .eq. 'y' .or. ans .eq. 'Y') then |
|
427 ! slice = .true. |
|
428 ! else |
|
429 ! slice = .false. |
|
430 ! end if |
|
431 ! |
|
432 slice = .false. |
|
433 ! |
|
434 ! print 104 |
|
435 ! 104 format (t10,'Addition <Y,N>?') |
|
436 ! read (5,105) ans |
|
437 ! 105 format (a1) |
|
438 ! if (ans .eq. 'y' .or. ans .eq. 'Y') then |
|
439 ! addition = .true. |
|
440 ! else |
|
441 ! addition = .false. |
|
442 ! end if |
|
443 addition = .false. |
|
444 ! |
|
445 ! call 2d multi-fractal generator |
|
446 print *,"Calling fractal generator for category ", iz |
|
447 call mpfm2d(X, maxlevel, sigma, H, addition, wrap, gradient, iz) |
|
448 ! |
|
449 print *, "fractal generated" |
|
450 xmin = 1e30 |
|
451 xmax = -1e30 |
|
452 xmean = 0.0 |
|
453 xss = 0.0 |
|
454 icnt = 0 |
|
455 ! |
|
456 ! obtain xmin and xmax for range normalization |
|
457 ! calculate initial raw statistics |
|
458 ! while you're at it |
|
459 do i = 0,N |
|
460 do j = 0,N |
|
461 if (X(i,j) .lt. xmin) xmin = X(i,j) |
|
462 if (X(i,j) .gt. xmax) xmax = X(i,j) |
|
463 ! xmean = xmean + X(i,j) |
|
464 ! xss = xss + X(i,j)*X(i,j) |
|
465 ! icnt = icnt+1 |
|
466 end do |
|
467 end do |
|
468 ! xmean = xmean / icnt |
|
469 ! print 200, xmean, xss, xmin, xmax |
|
470 !200 format (/t10,'Hfract Summary'/ & |
|
471 ! t10,'Mean = ',g12.6 & |
|
472 ! /t10,'Sum Sq. = ',g12.6 & |
|
473 ! /t10,'Minimum = ',g12.6 & |
|
474 ! /t10,'Maximum = ',g12.6) |
|
475 ! |
|
476 ! rescale and normalize fractal map between 0 and 32767, inclusive |
|
477 ! |
|
478 range1 = xmax - xmin |
|
479 do i = 0,N |
|
480 do j = 0,N |
|
481 iv = (i*(N+1))+j+1 |
|
482 P(iv,iz) = ((X(i,j) - xmin)/range1) * 32767 |
|
483 P(iv,ihabs+1) = i |
|
484 P(iv,ihabs+2) = j |
|
485 ! write(10,*) P(iv,iz) |
|
486 end do |
|
487 end do |
|
488 print*,'done normalizing now' |
|
489 ! |
|
490 !c write GRASS ascii map |
|
491 ! write(10,791), N |
|
492 !791 format ('north: ',i4) |
|
493 ! write(10,792) |
|
494 !792 format ('south: 0') |
|
495 ! write(10,793), N |
|
496 !793 format ('east: ',i4) |
|
497 ! write(10,794) |
|
498 !794 format ('west: 0') |
|
499 ! write(10,795), N+1 |
|
500 !795 format ('rows: ',i4) |
|
501 ! write(10,796), N+1 |
|
502 !796 format ('cols: ',i4) |
|
503 ! do i = 0,N |
|
504 ! write(10,*) (P((i*(N+1))+j+1,iz), j=0,N) |
|
505 ! end do |
|
506 ! |
|
507 ! end of iz loop over all landscape categories in synthetic map |
|
508 22 continue |
|
509 end do |
|
510 ! ********************************************************************** |
|
511 ! |
|
512 ! ihabs fractal maps now installed in P matrix: |
|
513 ! |
|
514 ! ihabs cols + 3 |
|
515 ! |---------------| |
|
516 ! | | |
|
517 ! | | |
|
518 ! | empty cells | |
|
519 ! mtresolved =>| | |
|
520 ! | | |
|
521 ! mtsplit =>|---------------| |
|
522 ! | | |
|
523 ! | | |
|
524 ! | | |
|
525 ! | | |
|
526 ! | singly- | |
|
527 ! | painted | |
|
528 ! | | |
|
529 ! | | |
|
530 ! | | |
|
531 ! | | |
|
532 ! | | |
|
533 ! itiedsplit =>|---------------| |
|
534 ! | | |
|
535 ! | | |
|
536 ! | ties | |
|
537 ! m =>| | |
|
538 ! | | |
|
539 ! numcells =>|_______________| |
|
540 ! |
|
541 ! ihabs+1 column carries x-coord of cell |
|
542 ! ihabs+2 column carries y-coord of cell |
|
543 ! ihabs+3 column carries number of ties for this cell |
|
544 ! |
|
545 ! |
|
546 ! initialize II array with background value for re-use as output map |
|
547 II(:,:) = 32767 |
|
548 ! |
|
549 ! set paint thresholds |
|
550 print *,"Calculating paint thresholds" |
|
551 do i = 0,ihabs ! over all map categories |
|
552 ! |
|
553 !c don't bother if no cells requested of this cat |
|
554 if (prob(i) .ne. 0.0) then |
|
555 ! |
|
556 ! istart = 0 |
|
557 ! print *, "numcells is ", numcells |
|
558 ! print *, "istart is ", istart |
|
559 ! print *, "i is ", i |
|
560 ! Call HeapSort(P,numcells,istart,i) ! sort P matrix by ith column |
|
561 ! |
|
562 ! |
|
563 ! next category to be painted goes to single vector for sorting |
|
564 idispvec(0) = 0 |
|
565 do l = 1, numcells |
|
566 idispvec(l) = P(l,i) |
|
567 end do |
|
568 ! |
|
569 ! call special vector version of Heapsort |
|
570 istart = 1 |
|
571 Call HS2(idispvec,numcells+1,istart) |
|
572 ! |
|
573 !c diagnostic print |
|
574 ! do k=16636,16646 |
|
575 ! print *, i, k, idispvec(k) |
|
576 ! end do |
|
577 ! |
|
578 ! set paint threshold for this map category |
|
579 hslice(i) = idispvec(numcells - (int((prob(i) * numcells)+.499))) |
|
580 ! for absolute constraint maps |
|
581 if (hslice(i) .eq. 32767 .and. constraint(i)) hslice(i) = 32766 |
|
582 print *, "trying to assign ", int((prob(i) * numcells)+.499), & |
|
583 " cells for category", i |
|
584 jsum = jsum + int((prob(i) * numcells)+.499) |
|
585 end if |
|
586 33 end do ! over all map categories |
|
587 ! |
|
588 print*, 'trying to assign a total of', jsum,' cells of', numcells, & |
|
589 ' total cells' |
|
590 !c |
|
591 !c set paint threshold for this map category |
|
592 ! hslice(i) = P(numcells - (int((prob(i) * numcells)+.499)),i) |
|
593 ! print *, "trying to assign ", int((prob(i) * numcells)+.499), & |
|
594 ! " cells for category", i |
|
595 ! end do ! over all map categories |
|
596 ! |
|
597 ! dump hslice array for checking |
|
598 print *, 'horizontal slices at these elevations:' |
|
599 do i = 0,ihabs |
|
600 print *, i, hslice(i) |
|
601 end do |
|
602 ! |
|
603 ! write XPM file of fractal probability terrain aspect |
|
604 if (probout) then |
|
605 do iz = 0, ihabs |
|
606 if (.not. (constraint(iz) .or. prob(iz) .eq. 0.0 .or. & |
|
607 iz .eq. igreatest)) then |
|
608 Call XPMRelief(P,iz,hslice(iz)) |
|
609 end if |
|
610 335 end do |
|
611 end if |
|
612 ! |
|
613 ! |
|
614 numpainted = 0 |
|
615 numtied = 0 |
|
616 numblank = 0 |
|
617 ! |
|
618 print *,"Beginning tie resolution process" |
|
619 do l = 1, numcells ! start at top of P matrix |
|
620 idummy = 0 |
|
621 do j = 0, ihabs ! over all map categories |
|
622 ! probably don't need to do this if prob(j) .eq. 0.0 |
|
623 if (.not. (prob(j) .eq. 0.0 .or. j .eq. igreatest)) then |
|
624 ! |
|
625 ! if this cell is painted this category |
|
626 if ( P(l,j) .gt. hslice(j) ) then |
|
627 ! if ( P(l,j) .ge. hslice(j) ) then |
|
628 ! |
|
629 ! count the multiple assignments |
|
630 P(l,ihabs+3) = P(l,ihabs+3) + 1 |
|
631 ! |
|
632 ! find the winner category for this cell |
|
633 ! winner category has the highest probability surface |
|
634 if ( P(l,j) .gt. idummy ) then |
|
635 ihighest = j |
|
636 idummy = P(l,ihighest) |
|
637 end if |
|
638 ! |
|
639 end if ! this cell painted this category |
|
640 ! |
|
641 end if |
|
642 end do ! over all map categories |
|
643 ! print *, "highest for cell ", l," was category", ihighest |
|
644 ! |
|
645 ! keep count of each type of cell for displacement into |
|
646 ! sorted array later |
|
647 ! if singly-painted or currently tied, paint outmap with ihighest |
|
648 if ( P(l,ihabs+3) .eq. 1) then |
|
649 ! paint this cell with the winning category |
|
650 ! print *, "painting uncontested cell", P(l,ihabs+1), & |
|
651 ! P(l,ihabs+2)," with cat", ihighest |
|
652 ! |
|
653 numpainted = numpainted + 1 |
|
654 II( P(l,ihabs+1),P(l,ihabs+2) ) = ihighest |
|
655 else if ( P(l,ihabs+3) .eq. 0) then |
|
656 numblank = numblank + 1 |
|
657 else |
|
658 numtied = numtied + 1 |
|
659 ! print *, "painting tied cell", P(l,ihabs+1),P(l,ihabs+2), & |
|
660 ! " with winning cat", ihighest |
|
661 II( P(l,ihabs+1),P(l,ihabs+2) ) = ihighest |
|
662 end if |
|
663 ! |
|
664 end do |
|
665 ! |
|
666 ! at this point, II array has all winners and uncontested |
|
667 ! cells painted |
|
668 ! |
|
669 ! sort P array into 3 sections by number of assignments |
|
670 ! ascending sort results in empty, singly-painted, and |
|
671 ! ties sections |
|
672 print *,"Preparing to sort ties for resolution" |
|
673 ! |
|
674 ! sort by ties column |
|
675 istart = 1 |
|
676 Call HeapSort(P,numcells,istart,ihabs+3) |
|
677 ! print *, "got back from HeapSort" |
|
678 ! |
|
679 ! calculate pointers to section breaks |
|
680 mtsplit = numblank |
|
681 itiedsplit = numblank + numpainted |
|
682 ! |
|
683 mtresolved = mtsplit |
|
684 ! |
|
685 print *, "numblank is ", numblank |
|
686 print *, "numtied is ", numtied |
|
687 print *, "itiedsplit is ", itiedsplit |
|
688 print *, "numpainted is ", numpainted |
|
689 print *, "mtresolved and mtsplit are ", mtresolved |
|
690 |
|
691 if (numtied .gt. numblank) print*, 'WARNING: More tied cells', & |
|
692 ' than leftover blank cells to use for tie resolution!' |
|
693 ! |
|
694 ! sum multiwayness of ties and compare to number of empty cells |
|
695 do m = itiedsplit+1, numcells |
|
696 isum = isum + ( P(m,ihabs+3) - 1 ) |
|
697 end do |
|
698 print *, "sum of multiwayness of ties is ", isum |
|
699 if (isum .gt. numblank) print *, 'Warning: cells are tied more', & |
|
700 ' ways than there are blank cells to assign!' |
|
701 ! |
|
702 ! call dump_map(P) |
|
703 ! |
|
704 ! sort ties into priority order for resolution |
|
705 ! by calculating and sorting on summed elevation |
|
706 ! for all tied categories |
|
707 ! |
|
708 do m = itiedsplit+1, numcells ! over all tied |
|
709 kount = 0 |
|
710 idummy = 0 |
|
711 do k = 0, ihabs ! over all categories |
|
712 ! |
|
713 ! probably not necessary if prob(k) .eq. 0 |
|
714 if (.not. (prob(k) .eq. 0.0 .or. k .eq. igreatest)) then |
|
715 ! if this cell is painted this category |
|
716 if ( P(m,k) .gt. hslice(k) ) then |
|
717 ! if ( P(m,k) .ge. hslice(k) ) then |
|
718 kount = kount + 1 ! number of ties for this cell |
|
719 ties(kount) = k ! keep categories of ties |
|
720 ! |
|
721 end if |
|
722 end if |
|
723 933 end do |
|
724 ! |
|
725 zmean = 0.0 |
|
726 do k = 1, kount |
|
727 zmean = zmean + P(m,ties(k)) |
|
728 end do |
|
729 ! P(m,ihabs+3) = zmean/kount |
|
730 ! try weighting priority by simple sum, not average |
|
731 ! therefore, max sum = (ihabs-1) * 32767 |
|
732 ! rescale as proportion of 32767 for integer*4 array cell |
|
733 P(m,ihabs+3) = zmean/(ihabs-1) |
|
734 ! |
|
735 end do |
|
736 ! |
|
737 ! sort ties by mean priority |
|
738 Call HeapSort(P,numtied,itiedsplit+1,ihabs+3) |
|
739 ! |
|
740 ! save min and max tie priorities for gray levels in XPMTies |
|
741 maxgray = P(numcells,ihabs+3) |
|
742 print*, "maxgray is ", maxgray |
|
743 ! |
|
744 ! generate xpm picture of tied cells |
|
745 ! colored by priority gray shades |
|
746 ! Call XPMTies(P) |
|
747 ! |
|
748 ! |
|
749 ! start at bottom of list (highest priority ties) |
|
750 do m = numcells, itiedsplit+1, -1 ! over all tied |
|
751 kount = 0 |
|
752 idummy = 0 |
|
753 do k = 0, ihabs ! over all categories |
|
754 ! |
|
755 ! probably not necessary if prob(k) .eq. 0 |
|
756 if (.not. (prob(k) .eq. 0.0 .or. k .eq. igreatest)) then |
|
757 ! if this cell is painted this category |
|
758 if ( P(m,k) .gt. hslice(k) ) then |
|
759 ! if ( P(m,k) .ge. hslice(k) ) then |
|
760 kount = kount + 1 ! number of ties for this cell |
|
761 ties(kount) = k ! keep categories of ties |
|
762 ! save most likely category for this cell |
|
763 if ( P(m,k) .gt. idummy ) then |
|
764 ihighest = k |
|
765 idummy = P(m,ihighest) |
|
766 end if |
|
767 end if |
|
768 end if |
|
769 656 end do |
|
770 ! |
|
771 ! |
|
772 ! resolve all ties - |
|
773 ! over all tied categories that are not the winner |
|
774 ! (which has already been painted) |
|
775 ! we need an empty cell - |
|
776 ! search the empty cell list by this cat to find next |
|
777 ! most likely available empty cell |
|
778 ! |
|
779 ! empty cell list will be searched once for each tie category |
|
780 ! except the winning category |
|
781 ! |
|
782 ! print *, "deciding among ", kount, " ties for cell ", m |
|
783 do k = 1, kount ! over all ties for this cell |
|
784 if ( ties(k) .ne. ihighest ) then ! don't re-do the winner |
|
785 ! print *, ties(k), " was a tied category for cell ", m |
|
786 ! |
|
787 ibiggest = 0 |
|
788 idummy = 0 |
|
789 do i = 1,mtresolved |
|
790 if (P(i,ties(k)) .gt. idummy) then |
|
791 ibiggest = i |
|
792 idummy = P(i,ties(k)) |
|
793 end if |
|
794 end do |
|
795 ! |
|
796 Call Interchange(P,ibiggest,mtresolved) |
|
797 ! |
|
798 ! istart = 1 |
|
799 ! Call HeapSort(P,mtresolved,istart,ties(k)) |
|
800 ! |
|
801 ! now mtresolved points to highest |
|
802 ! next available empty cell of category ties(k) |
|
803 ! paint it! |
|
804 ! print *, "painting empty cell", mtresolved," at coords ", & |
|
805 ! P(mtresolved,ihabs+1), P(mtresolved,ihabs+2)," with prob ", & |
|
806 ! P(mtresolved,ties(k))," with cat", ties(k), & |
|
807 ! " for tie at cell", m |
|
808 II( P(mtresolved,ihabs+1),P(mtresolved,ihabs+2) ) = ties(k) |
|
809 ! burn one empty |
|
810 mtresolved = mtresolved - 1 |
|
811 end if |
|
812 ! |
|
813 992 end do ! over all ties for this cell |
|
814 ! print *, "resolving ties at cell ", m |
|
815 ! |
|
816 end do ! over all tied cells |
|
817 ! |
|
818 !********************************************************************** |
|
819 ! |
|
820 ! |
|
821 ! Call dump_map(P) |
|
822 ! |
|
823 ! generate xpm picture of tied cells |
|
824 ! colored by priority gray shades |
|
825 ! and used blank cells colored in yellow |
|
826 Call XPMTies(P) |
|
827 ! |
|
828 ! |
|
829 ! set background category and |
|
830 ! perform summary of realized map |
|
831 do i = 0,N |
|
832 do j = 0,N |
|
833 if (II(i,j) .eq. 32767) II(i,j) = igreatest |
|
834 kcnt(II(i,j)) = kcnt(II(i,j)) + 1 |
|
835 ! write(12,*) II(i,j) |
|
836 end do |
|
837 end do |
|
838 ! |
|
839 ! if requested, |
|
840 ! write GRASS ascii map |
|
841 if (grassout) then |
|
842 write(10,91), N |
|
843 91 format ('north: ',i4) |
|
844 write(10,92) |
|
845 92 format ('south: 0') |
|
846 write(10,93), N |
|
847 93 format ('east: ',i4) |
|
848 write(10,94) |
|
849 94 format ('west: 0') |
|
850 write(10,95), N+1 |
|
851 95 format ('rows: ',i4) |
|
852 write(10,96), N+1 |
|
853 96 format ('cols: ',i4) |
|
854 write(fmt,'("(",I0,"(i5,1x))")') j |
|
855 do i = 0,N |
|
856 ! write(10,47) (II(i,j), j=0,N) |
|
857 write(10,fmt) (II(j,i), j=0,N) |
|
858 ! write(10,47) (II(j,i), j=0,N) |
|
859 !47 format (<j>(i5,1x)) |
|
860 end do |
|
861 end if |
|
862 ! |
|
863 ! write final map to xpm file |
|
864 if (xpmout) Call XPMMap(N,II) |
|
865 ! |
|
866 print 299 |
|
867 299 format (/t10,'Actual Map Summary:'/t12,'Habitat',5x, & |
|
868 'P',3x,'Real Cells',3x,'Real P') |
|
869 ! pfract = 0.0 |
|
870 do i = 0, ihabs |
|
871 pfract = real(kcnt(i),kind(pfract)) / real(numcells,kind(pfract)) |
|
872 print 300, i, prob(i), kcnt(i), pfract |
|
873 300 format (t10, i6, 4x, f6.5, 2x, i6, 2x, f6.5) |
|
874 end do |
|
875 ! |
|
876 !********************************************************************** |
|
877 ! |
|
878 ! visualize the synthetic landscape |
|
879 if (vizout) then |
|
880 ! |
|
881 ! print final map to screen |
|
882 do i=0,N |
|
883 do j=0,N |
|
884 idispvec(j) = II(i,j) |
|
885 end do |
|
886 ! print 199, (idispvec(j), j=0,N) |
|
887 !199 format (<N+1>(i1, 1x)) |
|
888 end do |
|
889 ! |
|
890 ! print current state of P array |
|
891 ! Call print(P) |
|
892 ! |
|
893 ! visualize II array |
|
894 mapcolor = 'bobs.cmap\0' |
|
895 maptitle = 'Fractal Landscape Realizer \0' |
|
896 idtype = 4 |
|
897 idumm = 0 |
|
898 ipixid = -1 |
|
899 ! |
|
900 ! build display vector |
|
901 ! make small maps more visible by repeating cells and lines |
|
902 ! irepeat = maxN / N |
|
903 ! print *, irepeat |
|
904 k = 0 |
|
905 do i = 0,N |
|
906 ! kstart = k |
|
907 do j=0,N |
|
908 ! do kk = 1, irepeat |
|
909 ! idispvec(k)=II(i,j) |
|
910 idispvec(k)=II(j,i) |
|
911 k = k+1 |
|
912 ! end do ! repeating cells in this row |
|
913 end do ! all cells in this row |
|
914 ! kend = k-1 |
|
915 !c repeat lines |
|
916 ! if (irepeat .ge. 2) then |
|
917 ! do jj = 1,irepeat-1 |
|
918 ! do kk = kstart, kend |
|
919 ! idispvec(k) = idispvec(kk) |
|
920 ! k = k+1 |
|
921 ! end do ! dup as many times as necessary |
|
922 ! end do ! over all repeated lines |
|
923 ! end if |
|
924 ! |
|
925 end do |
|
926 ! |
|
927 !c adjust rows and columns by multiplier |
|
928 ! isize = (N+1) * irepeat |
|
929 ! jsize = (N+1) * irepeat |
|
930 ! print *, isize, jsize |
|
931 ! ierr=DrawPixmap(ipixid, isize, jsize, idtype, idispvec, |
|
932 ! ierr=DrawPixmap(ipixid,N+1,N+1,idtype,idispvec,idumm,mapcolor, & |
|
933 ! maptitle) |
|
934 ! print *,'ierr is', ierr |
|
935 ! print *, 'done visualizing' |
|
936 ! call sleep(500) |
|
937 ! |
|
938 end if ! if visualizing landscape |
|
939 ! |
|
940 stop |
|
941 ! end program MidPtFM2d |
|
942 end program realizer |