-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathmaphandler.py
2565 lines (2289 loc) · 142 KB
/
maphandler.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# -*- coding: utf-8 -*-
#------------------------------------------------------------------------------
#
# Name: l maphandler.py
# Purpose: This module includes all functions for the usage of maps.
#
# Author: Carola Paetzold, Sven Lautenbach, Michael Strauch
# Contact: [email protected]
#
# Helmholtz Centre for Environmental Research - UFZ
# Department Computational Landscape Ecology - CLE
# Permoserstrasse 15
# D-04318 Leipzig, Germany
# http://www.ufz.de
#
# Created: We Oct 29 2014
#
# Copyright: (c) Carola Paetzold / Sven Lautenbach / Michael Strauch 2017
#
# Licence: This program is free software:
# you can redistribute it and/or modify it under the terms
# of the GNU General Public License as published by the
# Free Software Foundation, either version 3 of the License,
# or (at your option) any later version. This program is
# distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty
# of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
# See the GNU General Public License for more details.
# You should have received a copy of the GNU General
# Public License along with this program.
# If not, see <http://www.gnu.org/licenses/>.
#
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
# Imports
#------------------------------------------------------------------------------
import os
import sys
import csv
import time
import random
import numpy as np
from math import pow
import requirements as req
from filehandler import WriteLogMsg
from filehandler import WriteMap
from filehandler import WriteCandidateList
import config as cfg
wrkDir = os.path.abspath('.')
# first_ind is True for the first call of the generate_parameter function
first_ind = True
# array for the start individual
start_individual = []
# array for original ASCII map
map = []
# array for the patch ID map
patchID_map = []
# array for the header information from ASCII map
header = []
# list with header information from ASCII map - one string per line
header_all = []
# array for land use transition constraints
trans_matrix = []
# array for total area constraints of each land use class
min_max_diff = []
# array for static land use types
static_elements = []
# array for non static land use types
nonstatic_elements = []
# array for the area percentage of the patches
map_proportion = []
# dictionary for possible land use options per element
possible_elements = {}
# array with segments of impossible candidates
impossible_cand = []
# Boolean variable for starting the special_termination
end_optimization = False
# for constraint tournament selection
# dictionary as memory for the analysis results of the candidate violation
violation_memory = {}
# dictionary for area percentage of the static land use categories
static_area = {}
# dictionary with patch indices and total area per land use class
start_land_cover = {}
# dictionary with information per land use class if extreme seeds were created before
extreme_seeds_dict = {}
#-------------------------------------------------------------------------------------
# Start the termination of the optimization algorithm
#-------------------------------------------------------------------------------------
def logical_termination():
"""If repair_mutation is selected and the last possible candidate is set as an individual
then the variable end = True and the termination can start after the determination of
the fitness values of the last population.
"""
# Boolean variable for starting the special_termination
global end_optimization
return end_optimization
#-------------------------------------------------------------------------------------
# Read an ASCII-map
#-------------------------------------------------------------------------------------
def read_ascii_map(file):
"""Read an ascii file.
Return the map as matrix and header data as array.
input data: ascii file
"""
msg = "Read ascii file ..."
WriteLogMsg(msg)
# read header information in the header array
header_all = open(file,'rb').readlines()[0:6]
header = []
for line in header_all:
line_split = line.split()
header.append(line_split[1])
# read the map in the map matrix
map = np.genfromtxt(file, dtype=int, skip_header=6, filling_values='-1')
print("map \n%s" %map)
return header, header_all, map
#-------------------------------------------------------------------------------------
# Determine the coordinates of neighboring cells of cell (col, row)
#-------------------------------------------------------------------------------------
def getNbh(col, row, ncols, nrows, four_neighbours):
"""Determine the neighboring cells of the cell (col,row) and
return the coordinates as arrays separated in nbhs_col and nbhs_row.
The combination of the elements gives the coordinates of the neighbouring cells.
input:
col and row are coordinates of the reviewed element
ncols, nrows are numbers of rows and columns in the map
four_neighbours if True than 4 neighboring cells are scanned else 8
"""
# assuming that a cell in the center has 8 neighbouring cells
if four_neighbours == 'False':
# cell is no edge cell
if col > 0 and row > 0 and row < nrows -1 and col < ncols -1:
nbhs_col = [x + col for x in[-1, -1, -1, 0, 0, 1, 1, 1]]
nbhs_row = [x + row for x in[-1, 0, 1, -1, 1, -1, 0, 1]]
# cell is a left edge element but no corner element
elif col == 0 and row > 0 and row < nrows -1:
nbhs_col= [x + col for x in[0, 1, 1, 0, 1]]
nbhs_row= [x + row for x in[-1, -1, 0, 1, 1]]
# cell is a right edge element but no corner element
elif col == ncols -1 and row > 0 and row < nrows -1:
nbhs_col= [x + col for x in[-1, -1, -1, 0, 0]]
nbhs_row= [x + row for x in[-1, 0, 1, -1, 1]]
# cell is an upper edge element but no corner element
elif row == 0 and col > 0 and col < ncols -1:
nbhs_col= [x + col for x in[-1, -1, 0, 1, 1 ]]
nbhs_row= [x + row for x in[ 0, 1, 1, 0, 1 ]]
# cell is a bottom edge element but no corner element
elif row == nrows -1 and col > 0 and col < ncols -1:
nbhs_col= [x + col for x in[-1, -1, 0, 1, 1 ]]
nbhs_row= [x + row for x in[ -1, 0, -1, -1, 0 ]]
# cell is in the left upper corner
elif col == 0 and row == 0:
nbhs_col= [x + col for x in[ 0, 1, 1]]
nbhs_row= [x + row for x in[ 1, 0, 1]]
# cell is in the left bottom corner
elif col == 0 and row == nrows -1:
nbhs_col= [x + col for x in[ 0, 1, 1]]
nbhs_row= [x + row for x in[ -1, 0, -1]]
# cell is in the right upper corner
elif col == ncols -1 and row == 0:
nbhs_col= [x + col for x in[ -1, -1, 0]]
nbhs_row= [x + row for x in[ 0, 1, 1]]
# cell is in the right bottom corner
else:
nbhs_col= [x + col for x in[ -1, -1, 0 ]]
nbhs_row= [x + row for x in[ -1, 0, -1]]
# assuming that a cell in the center has 4 neighbouring cells
elif four_neighbours == 'True':
# cell is no edge cell
if col > 0 and row > 0 and row < nrows -1 and col < ncols -1:
nbhs_col = [x + col for x in[-1, 0, 0, 1]]
nbhs_row = [x + row for x in[ 0, -1, 1, 0]]
# cell is a left edge element but no corner element
elif col == 0 and row > 0 and row < nrows -1:
nbhs_col= [x + col for x in[0, 1, 0]]
nbhs_row= [x + row for x in[-1, 0, 1]]
# cell is a right edge element but no corner element
elif col == ncols -1 and row > 0 and row < nrows -1:
nbhs_col= [x + col for x in[-1, 0, 0]]
nbhs_row= [x + row for x in[ 0, 1, -1]]
# cell is an upper edge element but no corner element
elif row == 0 and col > 0 and col < ncols -1:
nbhs_col= [x + col for x in[-1, 0, 1]]
nbhs_row= [x + row for x in[ 0, 1, 0]]
# cell is an bottom edge element but no corner element
elif row == nrows -1 and col > 0 and col < ncols -1:
nbhs_col= [x + col for x in[-1, 0, 1]]
nbhs_row= [x + row for x in[ 0, -1, 0]]
# cell is in the left upper corner
elif col == 0 and row == 0:
nbhs_col= [x + col for x in[ 0, 1]]
nbhs_row= [x + row for x in[ 1, 0]]
# cell is in the left bottom corner
elif col == 0 and row == nrows -1:
nbhs_col= [x + col for x in[ 0, 1]]
nbhs_row= [x + row for x in[ -1, 0]]
# cell is in the right upper corner
elif col == ncols -1 and row == 0:
nbhs_col= [x + col for x in[ -1, 0]]
nbhs_row= [x + row for x in[ 0, 1]]
# cell is in the right bottom corner
else:
nbhs_col= [x + col for x in[ -1, 0 ]]
nbhs_row= [x + row for x in[ 0, -1]]
else:
msg = "Error: ini input for four_neighbours is not correct. Please check."
WriteLogMsg(msg)
raise SystemError("Error: ini input for four_neighbours is not correct")
req.close_window
return [nbhs_row, nbhs_col]
#-------------------------------------------------------------------------------------
# Determination of patch elements
#-------------------------------------------------------------------------------------
def determine_patch_elements(row, col, map, patch_map, patch_ID, cls, four_neighbours):
"""This recursive function scans all patch elements
and returns the coordinates of these elements.
input:
col and row are coordinates of the parent element
map is the original ascii map
patch_map is a map with patch_IDs for each patch element
patch_ID is the ID of the new patch
cls is the land use index of the patch
four_neighbours if True than 4 neighboring cells are scanned else 8
"""
# determine coordinates of neighboring cells
new_nbhs_row, new_nbhs_col = getNbh(col, row, map.shape[1], map.shape[0], four_neighbours)
# stack for patch elements whose neighboring cells should be determined
nbhs_row = []
nbhs_col = []
for i in range(len(new_nbhs_row)):
# add new neighboring cells to nbhs_row/col if new cells belong to cls and are not jet marked as patch element
# the cell is no patch element if it has another land use id
if map[new_nbhs_row[i], new_nbhs_col[i]] == cls and patch_map[new_nbhs_row[i], new_nbhs_col[i]] == 0:
nbhs_row.append(new_nbhs_row[i])
nbhs_col.append(new_nbhs_col[i])
while len(nbhs_row) > 0:
# cells could be double in nbhs_row/col
if patch_map[nbhs_row[0], nbhs_col[0]] == 0:
# mark all patch elements in patch_map with patch_ID
patch_map[nbhs_row[0], nbhs_col[0]] = patch_ID
# get coordinates of neighboring cells of this cell
new_nbhs_row, new_nbhs_col = getNbh(nbhs_col[0], nbhs_row[0], map.shape[1], map.shape[0], four_neighbours)
for i in range(len(new_nbhs_row)):
# add new neighboring cells to nbhs_row/col if new cells belong to cls and are not jet marked as patch element
if map[new_nbhs_row[i], new_nbhs_col[i]] == cls and patch_map[new_nbhs_row[i], new_nbhs_col[i]] == 0:
nbhs_row.append(new_nbhs_row[i])
nbhs_col.append(new_nbhs_col[i])
# delete this checked neighboring cell of the array
del nbhs_row[0]
del nbhs_col[0]
return patch_map
#-------------------------------------------------------------------------------------
# Determine patch elements of a patch ID map and check equality of land use index
#-------------------------------------------------------------------------------------
def determine_IDmap_patch_elements(row, col, patch_map, map, neighbors, cls, landuse, error_dic, four_neighbours):
"""This recursive function scans all patch elements of the patch ID map,
check if all patch elements have the same land use index
and return the coordinates of the patch elements.
input:
col and row are coordinates of the parent element
patch_map is the given patch ID map
map is the original ascii map
neighbors is a matrix for marking the scanned cells
cls is the patch index
landuse is the index of the first scanned patch element
four_neighbours if True than 4 neighboring cells are scanned else 8
"""
# determine coordinates of neighboring cells
new_nbhs_row, new_nbhs_col = getNbh(col, row, map.shape[1], map.shape[0], four_neighbours)
# stack for patch elements whose neighboring cells should be determined
nbhs_row = []
nbhs_col = []
for i in range(len(new_nbhs_row)):
# add new neighboring cells to nbhs_row/col if new cells belong to cls in patch_map and are not jet marked as scanned
if patch_map[new_nbhs_row[i], new_nbhs_col[i]] == cls and neighbors[new_nbhs_row[i], new_nbhs_col[i]] == 0:
nbhs_row.append(new_nbhs_row[i])
nbhs_col.append(new_nbhs_col[i])
# print ("nbhs_row, nbhs_col von (%s,%s): %s, %s" %(row, col, nbhs_row, nbhs_col))
while len(nbhs_row) > 0:
# if cell was not scanned before
if neighbors[nbhs_row[0], nbhs_col[0]] == 0:
# mark this cell in neighbors with True (scanned)
neighbors[nbhs_row[0], nbhs_col[0]] = 1
# and check if land use type for all patch elements is equal
if map[nbhs_row[0], nbhs_col[0]] != landuse:
error_dic.update({"(%s,%s)" %(nbhs_row[0], nbhs_col[0]) : "more than one land use index for one patch" })
# determine coordinates of neighboring cells from this cell
new_nbhs_row, new_nbhs_col = getNbh(nbhs_col[0], nbhs_row[0], map.shape[1], map.shape[0], four_neighbours)
for i in range(len(new_nbhs_row)):
# add new neighboring cells to nbhs_row/col if new cells belong to cls in patch_map and not jet marked as scanned
if patch_map[new_nbhs_row[i], new_nbhs_col[i]] == cls and neighbors[new_nbhs_row[i], new_nbhs_col[i]] == 0:
nbhs_row.append(new_nbhs_row[i])
nbhs_col.append(new_nbhs_col[i])
# delete this checked neighboring cell of the array
del nbhs_row[0]
del nbhs_col[0]
return neighbors, error_dic
#-------------------------------------------------------------------------------------
# Cluster the cells of the map into patches
#-------------------------------------------------------------------------------------
def create_patch_ID_map(map, NODATA_value, static_elements, four_neighbours):
"""This function clusters the cells of the original map into patches
and returns a patch ID map as a 2 dimensional array and the start individual as vector.
input:
map is the original ascii map
NODATA_value is the NODATA_value of the original map
static_elements are the land use indices excluded from the optimization
four_neighbours if True than 4 neighboring cells are scanned else 8
"""
msg = "Cluster cells of the original map into patches ..."
WriteLogMsg(msg)
begin = time.time()
patches= np.zeros([map.shape[0], map.shape[1]], int)
ids = 0
NoData = int(NODATA_value)
genom = []
# loop over all cells
for row in range(0, map.shape[0]):
if (row + 0.0) % (round(map.shape[0] / 10.0)) == 0 and (time.time() - begin) > 2 :
progress = ((row+0.0) / map.shape[0]) * 100
WriteLogMsg("progress for scanning the rows in percent: %s" %progress)
for col in range(0, map.shape[1]):
# patchID = 0 used for static_elements
# map element was not scanned before as patch element and is not a static element or the NODATA_value
if patches[row,col]==0 and static_elements.count(map[row, col])==0 and map[row,col]!=NoData:
cls = map[row, col]
# increment scanned patch ID
ids += 1
# marke this cell as scanned patch element
patches[row, col] = ids
determine_patch_elements(row,col, map, patches, ids, cls, four_neighbours)
# add the map cell value to the individual vector
genom.append(cls)
end = time.time()
WriteLogMsg("Runtime for the generation of the patch_ID_map: %d seconds." %(end-begin))
return patches, genom
#-------------------------------------------------------------------------------------
# Read the patch ID map and create the start individual
#-------------------------------------------------------------------------------------
def read_patch_ID_map(file, map, NODATA_value, static_elements, four_neighbours):
"""This function reads a given patch ID map, checks its plausibility
and returns the patch ID map as a 2 dimensional array and the start individual as vector.
input:
file with the patch ID map (ascii format)
map is the original ascii map
NODATA_value is the NODATA_value of the original map
static_elements are the land use indices excluded from the optimization
four_neighbours if True than 4 neighboring cells are scanned else 8
"""
msg = "Read the patch ID map ..."
WriteLogMsg(msg)
# transform map into a matrix
try:
patches = np.genfromtxt(file, dtype=int, skip_header=6, filling_values='-1')
except ValueError as e:
msg = "Error: Number of values in one row or column is not correct. Please check.\n %s" %e
WriteLogMsg(msg)
raise SystemError("Error: Number of values in one row or column is not correct. Please check.\n %s" %e)
req.close_window
# check if number of rows and columns are equal to the original map
if (map.shape[0] != patches.shape[0]) or (map.shape[1] != patches.shape[1]):
msg = "Error: Number of rows or columns of the original map (rows: %s, columns: %s) and the patch ID map (rows: %s, columns: %s) are not equal. Please check." %(map.shape[0],map.shape[1],patches.shape[0],patches.shape[1])
WriteLogMsg(msg)
raise SystemError("Error: Number of rows or columns of the original map and the patch ID map are not equal.")
req.close_window
# check that land use indices of the patch elements are equal
# and that the land use index is not a static element
# if all checks are okay then add the land use index to the genome of the start individual
max_patchID = patches.max()
genom = np.zeros(max_patchID, int)
help_map = np.zeros([map.shape[0], map.shape[1]], bool) # default is 0/False
error_dic = {}
for row in range(0, help_map.shape[0]):
for col in range(0, help_map.shape[1]):
# map element was not scanned before
if help_map[row,col]==0:
# first case: patch element = 0
# check if the index in the original map is a static element or NODATA_value
if patches[row,col]==0:
if static_elements.count(map[row, col])==0 and map[row,col]!=int(NODATA_value):
error_dic.update({"(%s,%s)" %(row, col) : "zero but no static element and not a NODATA_value" })
# patch element != 0
else:
# check that cell is not a static element or NODATA_value
if static_elements.count(map[row, col])!=0 or map[row,col]==int(NODATA_value):
error_dic.update({"(%s,%s)" %(row, col) : "non-zero but static element or NODATA_value" })
# if file_HRU == None then check that patch ID was not registered in the genome before
elif cfg.mapConfig.file_HRU == 'None' and genom[patches[row,col]-1] != 0:
error_dic.update({"(%s,%s)" %(row, col) : "non-zero but patch ID is used for more than one patch" })
# if file_HRU != None then check if land use of the two patches are equal
elif cfg.mapConfig.file_HRU != 'None' and genom[patches[row,col]-1] != 0 and genom[patches[row,col]-1] != map[row,col]:
error_dic.update({"(%s,%s)" %(row, col) : "non-zero but patch ID is used for more than one land use" })
# plausibility checks are okay, then
else:
# add the land use index to the genome
if genom[patches[row,col]-1] == 0:
genom[patches[row,col]-1] = map[row,col]
# scan all patch elements
# check if land use is equal
# and mark them in the help_map
determine_IDmap_patch_elements(row, col, patches, map, help_map, patches[row,col], genom[patches[row,col]-1], error_dic, four_neighbours)
# mark the scanned cell in help_map
# so the cell is not scanned again by the determine_IDmap_patch_elements
help_map[row,col] = 1
# if at least one plausibility check was not succesfull
# print the error messages and stop the program
if len(error_dic) != 0:
msg = "Error: Some cells of the patch ID map break the plausibility rules. Please check."
WriteLogMsg(msg)
for key in error_dic:
msg = "%s %s" %(key,error_dic[key])
WriteLogMsg(msg)
raise SystemError("Error: Some cells of the patch ID map break the plausibility rules.")
req.close_window
return patches, genom.tolist()
#-------------------------------------------------------------------------------------
# Check that original map and transition matrix are correct
#-------------------------------------------------------------------------------------
def check_matrices(map, trans_matrix):
"""This function checks the trans_matrix for the following conditions to be true:
- all diagonal elements are ones
- all elements except the first column and row are only ones and zeros
- first row and column elements are equal
- all indices of the map are included in the transition matrix
This function checks the map for the following conditions to be true::
- no elements are set by the default -1
If one condition is not fulfill the program ends with an error message.
input:
map is the original ascii map
trans_matrix holds the land use transition constraints
"""
if len(trans_matrix) != 0:
# check if all diagonal elements are ones
for row in range(1,trans_matrix.shape[0]):
if trans_matrix[row][row] != 1:
msg = "Error: Not all diagonal elements of the transition matrix are ones. Please check."
WriteLogMsg(msg)
raise SystemError("Error: Not all diagonal elements of the transition matrix are ones.")
req.close_window
# check if all elements except the first column and row are only ones and zeros
if not (len(np.unique(trans_matrix[1:,1:]))==2 and 1 in np.unique(trans_matrix[1:,1:]) and 0 in np.unique(trans_matrix[1:,1:])):
msg = "Error: Not all elements of the transition matrix (excluding first column and row) are ones and zeros. Please check."
WriteLogMsg(msg)
raise SystemError("Error: Not all elements of the transition matrix (excluding first column and row) are ones and zeros.")
req.close_window
# check if first row and column elements are equal
for row in range(1,trans_matrix.shape[0]):
if trans_matrix[row][0] != trans_matrix[0][row]:
msg = "Error: Elements of first row and column are not equal. Please check."
WriteLogMsg(msg)
raise SystemError("Error: Elements of first row and column are not equal.")
req.close_window
# no elements are set by the default -1
if (-1 in np.unique(trans_matrix)):
msg = "Error: Elements of the transition matrix are -1 (default value). Please check."
WriteLogMsg(msg)
raise SystemError("Error: Elements of the transition matrix are -1 (default value).")
req.close_window
# check if all indices of the map are included in the transition matrix
for element in np.unique(map):
if element not in np.unique(trans_matrix):
msg = "Error: Map element %s is not included in the transition matrix. Please check." % element
WriteLogMsg(msg)
raise SystemError("Error: Map element %s is not included in the transition matrix." %element)
req.close_window
# no elements are set by the default -1
if (-1 in np.unique(map)):
msg = "Error: Elements of the original map are -1 (default value). Please check."
WriteLogMsg(msg)
raise SystemError("Error: Elements of the original map are -1 (default value).")
req.close_window
# no elements of the original map are zero (no data < -1 and land use id > 0)
# 0 reserved for patch id map (patches which are excluded from the optimization)
if (0 in np.unique(map)):
msg = "Error: Elements of the original map are 0. No data should be < -1 and land use classes > 0. Please check."
WriteLogMsg(msg)
raise SystemError("Error: Elements of the original map are 0. No data should be < -1 and land use classes > 0. Please check.")
req.close_window
#-------------------------------------------------------------------------------------
# Determine all classes which are excluded from optimization
#-------------------------------------------------------------------------------------
def determine_static_classes(trans_matrix, max_range):
"""This function determines all classes which are excluded from optimization (static elements)
and returns arrays with the indices of static and non static elements.
input:
trans_matrix holding the land use transition constraints
max_range is the maximum number of possible land use options
"""
msg = "Determine static elements of the transition matrix ..."
WriteLogMsg(msg)
# identify all classes where column and row elements are zero (apart from diagonal elements)
# that means that all patches of these classes cannot be converted
static_elements = []
nonstatic_elements = []
# filter columns which fulfill the condition
ones = 0
# row-wise check for ones
for row in range(1,trans_matrix.shape[0]):
for col in range(1,trans_matrix.shape[1]):
if trans_matrix[row][col] == 1:
ones += 1
# mark the candidate as static or non static element (row is checked)
# if ones for row = 1 or max_range < land use index of trans_matrix
if ones==1 or trans_matrix[row][0] > max_range:
static_elements.append(trans_matrix[row][0])
else:
nonstatic_elements.append(trans_matrix[row][0])
ones = 0
# column-wise check for ones
ones = 0
index = 0
if len(static_elements) != 0:
for col in range(1,trans_matrix.shape[1]):
if index < len(static_elements) and static_elements[index] <= max_range:
if trans_matrix[0][col] == static_elements[index]:
for row in range(1,trans_matrix.shape[0]):
if trans_matrix[row][col] == 1:
ones += 1
if ones!=1:
# index remains as it is for the next round
# because of length reduction from static_element
nonstatic_elements.append(static_elements[index])
del static_elements[index]
else:
index += 1
ones = 0
if len(static_elements) == 0:
break
return static_elements, nonstatic_elements
#------------------------------------------------------------------------------
# Read the start individual and area proportion of each HRU from the HRU input file
#------------------------------------------------------------------------------
def read_HRUs(file_HRU):
"""The function counts the HRUs listed in the HRU input file
and returns the result.
input:
filename of the HRU list
"""
msg = "Read HRU file ..."
WriteLogMsg(msg)
# array for area percentage and name of the HRUs
global map_proportion
# set default value for areas excluded from optimization
map_proportion.append(float(0))
genom = []
reader = csv.reader(open(os.path.join(wrkDir, 'input', file_HRU), "r"))
i = 0
sum_area = 0
for row in reader:
if i != 0:
map_proportion.append(float(row[2]))
sum_area += float(row[2])
genom.append(int(row[0]))
i += 1
if round(sum_area,1) <= 100.0:
map_proportion[0] = 100 - sum_area
else:
msg = "Error: The input data for Area_rel are %s instead of 100 percent. Please check." %sum_area
WriteLogMsg(msg)
raise SystemError("Error: The input data for Area_rel are %s instead of 100 percent. Please check." %sum_area)
req.close_window
# save start individual in global variable
global start_individual
start_individual = genom
return genom
#-------------------------------------------------------------------------------------
# Generate the genome of the start individual
#-------------------------------------------------------------------------------------
def generate_genom(max_range, file_HRU, map_file, trans_file, patchIDmap_file, four_neighbours):
"""The function generates and returns the start individual and the non static land use indices
based on an HRU file or ASCII map.
It is called from generate_parameter in optiAlgorithm.py
and calls functions of maphandler.py for the individual generation.
input data:
max_range is the maximum number of possible land use options for the individual
file_HRU file with the original HRU list
map_file file with the original ascii map
trans_file file with the land use transition constraints
patchIDmap_file file with the patch ID map (ascii format)
four_neighbours if True than 4 neighboring cells are scanned else 8
"""
global map
global header
global header_all
global trans_matrix
global static_elements
global nonstatic_elements
global patchID_map
global map_proportion
# dictionary for possible land use options per element
global possible_elements
# HRU and given patch ID map need an original map
if file_HRU != 'None' and (map_file != 'None' or patchIDmap_file != 'None'):
if map_file == 'None' or patchIDmap_file == 'None':
msg = "Error: If you want to use maps with HRUs than you need an original map and a patch ID map. Please check."
WriteLogMsg(msg)
raise SystemError("Error: If you want to use maps with HRUs than you need an original map and a patch ID map. Please check.")
req.close_window
# if HRU file is given
if file_HRU != 'None':
msg = "Generate start individual based on an HRU file ..."
WriteLogMsg(msg)
genom = read_HRUs(file_HRU)
# if ASCII map is given
if map_file != 'None':
if file_HRU == 'None':
msg = "Generate start individual based on an ASCII map..."
WriteLogMsg(msg)
# read the original ascii map
header, header_all, map = read_ascii_map(os.path.join(wrkDir, 'input', map_file))
# no ASCII original map or HRU list is given
if map_file == 'None' and file_HRU == 'None':
msg = "Error: No original ASCII map or HRU list is given by the config.ini. Please check."
WriteLogMsg(msg)
raise SystemError("Error: No original ASCII map or HRU list is given by the config.ini.")
req.close_window
# if a transition matrix is given
if trans_file != 'None':
# read the transition matrix for land use change
msg = "Read transition matrix..."
WriteLogMsg(msg)
# write information into a matrix
trans_matrix = np.genfromtxt(os.path.join(wrkDir, 'input',trans_file), dtype=int, filling_values='-1')
# check if max_range is plausible
if max_range > trans_matrix.max():
msg = "Error: Input of max_range is bigger than the maximum element of the transition matrix. Please check."
WriteLogMsg(msg)
raise SystemError("Error: Input of max_range is bigger than the maximum element of the transition matrix.")
req.close_window
if map_file != 'None':
if max_range < map.max():
msg = "Error: Input of max_range is smaller than the maximum element (%s) in the original ASCII map. Please check." %map.max()
WriteLogMsg(msg)
raise SystemError("Error: Input of max_range is smaller than the maximum element in the original ASCII map.")
req.close_window
# check format requirements
check_matrices(map, trans_matrix)
# determine static land use elements
static_elements, nonstatic_elements = determine_static_classes(trans_matrix, max_range)
# possible land use options for each land use class according to transition matrix
if len(possible_elements) == 0:
for i in range(1,cfg.modelConfig.max_range+1):
if i not in static_elements:
# determine indices of transition matrix where in row with first element = i are ones
indices = np.asarray(np.nonzero(trans_matrix[np.nonzero(np.unique(trans_matrix[:,:1]) == i)[0][0],:])[0])
# delete the index for 1 as land use name
if 0 in indices:
indices = np.delete(indices,0)
# determine the land use classes of the indices
if len(indices) !=0 and trans_matrix[0][indices[0]] <= cfg.modelConfig.max_range:
land_use_classes = np.array([trans_matrix[0][indices[0]]])
for k in range(1,len(indices)):
if trans_matrix[0][indices[k]] <= cfg.modelConfig.max_range:
land_use_classes = np.append(land_use_classes,trans_matrix[0][indices[k]])
# save the possible land use classes
possible_elements.update({ i : land_use_classes})
WriteLogMsg("possible_elements per land use: ")
for key in possible_elements:
WriteLogMsg("land use, elements: %s,%s" %(key,possible_elements[key]))
# if no transition matrix is given -> land use change is completely random
else:
# check format requirements
if map_file != 'None':
check_matrices(map, trans_matrix)
# determine non_static elements ( 1 <= x <= max_range)
for i in range(1,max_range+1):
nonstatic_elements.append(i)
# header[5] is the NODATA_value of the original ascii map
# read patchID_map if it is available
if patchIDmap_file != 'None':
# read the patch ID map, check plausibility and create the start individual
if file_HRU == 'None':
msg = "Create the start individual based on the ID map ..."
WriteLogMsg(msg)
patchID_map, genom = read_patch_ID_map(os.path.join(wrkDir, 'input', patchIDmap_file), map, header[5], static_elements, four_neighbours)
# the patch ID map cannot be created based on the original map for HRUs (one HRU can consist of more than one patch)
elif file_HRU == 'None':
# create the patch ID map based on the original map and the start individual
patchID_map, genom = create_patch_ID_map(map, header[5], static_elements, four_neighbours)
# log the patch ID map
if file_HRU == 'None' or (file_HRU != 'None' and patchIDmap_file != 'None'):
WriteMap(header_all, patchID_map)
# analyze the proportional area of the patches in %
# for HRUs the map_proportion is obtained in read_HRUs
if file_HRU == 'None':
unique, counts = np.unique(patchID_map, return_counts = True)
number_cells = patchID_map.shape[0]*patchID_map.shape[1]
if 0 not in unique:
# no cells are excluded from optimization
map_proportion.append(float(0))
for item in counts:
map_proportion.append(float(item)/float(number_cells)*100)
# check if sum of area proportions is 100%
sum = 0
for i in map_proportion:
sum = sum + i
if round(sum) != 100:
msg = "Error: The sum of patch/HRU areas is %s instead of 100 percent. Please check." %sum
WriteLogMsg(msg)
raise SystemError("Error: The sum of patch/HRU areas is %s instead of 100 percent." %sum)
req.close_window
# check all indices of the map are included in the transition matrix
if trans_file != 'None':
for element in np.unique(genom):
if element not in np.unique(trans_matrix):
msg = "Error: Land use class %s is not included in the transition matrix. Please check." % element
WriteLogMsg(msg)
raise SystemError("Error: Land use class %s is not included in the transition matrix." %element)
req.close_window
# log the start individual and the non static land use indices
if len(genom) > 100:
msg = "start individual of the input data with length %d" %len(genom)
WriteLogMsg(msg)
else:
msg = "start individual of the input data with length %d: %r" %(len(genom),genom)
WriteLogMsg(msg)
msg = "DiscreteBounder values: %r" %nonstatic_elements
WriteLogMsg(msg)
# save start individual in global variable
global start_individual
start_individual = genom
WriteLogMsg("map proportion of the patches: %r" %map_proportion)
return genom, nonstatic_elements
#------------------------------------------------------------------------------
# Check if new_cand is element of the tabu memory
#------------------------------------------------------------------------------
def check_impossible_candidates(new_cand):
"""Return True if new_cand is element of the tabu memory
else False.
input data:
new_cand is a part of a new candidate
"""
# array with segments of impossible candidates (tabu memory)
global impossible_cand
# array for the start individual
global start_individual
return_value = False
copy_new_cand = new_cand[:]
# check complete candidate
if copy_new_cand in impossible_cand:
return_value = True
# check if first part is in impossible_cand:
else:
for item in impossible_cand:
if len(item) < len(copy_new_cand):
for index in range(len(item)):
# item not first part of new_cand
if item[index] != copy_new_cand[index]:
break
# item is first part of new_cand -> new_cand added to tabu memory
elif index == (len(item)-1) and item[index] == copy_new_cand[index]:
return_value = True
return return_value
#------------------------------------------------------------------------------
# Check if created individual is plausible
#------------------------------------------------------------------------------
def individual_filter(new_cand):
"""Check if the created individual violates the constraints.
Return True if the individual is feasible else False.
The function is structured for changeable objects for comparison but at
the moment only the original start individual is used for the comparison
and the decision if the new candidate is feasible or not.
input data:
new_cand is new candidate for checking the feasibility
"""
return_value = True
# array for land use transition constraints
global trans_matrix
# array for the area percentages of the patches
global map_proportion
# array for total area constraints
global min_max_diff
# array for the start individual
global start_individual
# array for static land use types
global static_elements
compare_individual = start_individual
# for each element of the individual, check if transition is allowed
if cfg.mapConfig.file_transformation != 'None':
for i in range(0,len(new_cand)):
if trans_matrix[np.nonzero(np.unique(trans_matrix[:,:1]) == compare_individual[i])[0][0]][np.nonzero(trans_matrix[0] == new_cand[i])[0][0]] != 1:
return_value = False
break
# check if total area constraints are satisfied
if cfg.mapConfig.file_difference != 'None':
for i in np.unique(new_cand):
sum_new = 0
for m in np.nonzero(new_cand==i)[0]:
# index m+1 because of the exclusion of patch ID 0 from optimization
sum_new += map_proportion[m+1]
if not (min_max_diff[1][np.nonzero(min_max_diff[0] == i)[0][0]] <= sum_new <= min_max_diff[2][np.nonzero(min_max_diff[0] == i)[0][0]]):
return_value = False
break
# check if area rules for land use types in new_cand do not apply -> min should be zero
for i in range(1,cfg.modelConfig.max_range+1):
if (i not in static_elements) and (i not in np.unique(new_cand)):
if min_max_diff[1][np.nonzero(min_max_diff[0] == i)[0][0]] > 0:
return_value = False
break
return return_value
#------------------------------------------------------------------------------
# Determine possible land use classes for next new_cand element
#------------------------------------------------------------------------------
def determine_possible_elements(possible_land_use, new_cand):
"""This function is used during repair mutation and returns the
possible land use classes for the next new_cand element taking
into account the tabu memory.
input data:
possible_land_use is an array with possible land use classes
derived from the transition matrix
new_cand is a part of a new candidate
"""
# array with segments of impossible candidates (tabu memory)
global impossible_cand
copy_poss_el = np.copy(possible_land_use)
copy_new_cand = new_cand[:]
# create all new candidate parts which can be created with one of the possible_elements and the new_cand
# and check if the new one is excluded by the tabu memory
for element in possible_land_use:
if len(copy_new_cand) == len(new_cand):
copy_new_cand.append(element)
else:
copy_new_cand[len(copy_new_cand)-1] = element
if copy_new_cand in impossible_cand:
copy_poss_el = np.delete(copy_poss_el,np.nonzero(copy_poss_el == copy_new_cand[len(copy_new_cand)-1])[0][0])
return copy_poss_el
#------------------------------------------------------------------------------
# Determine possible land use classes for next new_cand element
#------------------------------------------------------------------------------
def determine_possible_landuse(possible_land_use, index, new_cand):
"""This function is used during extreme seed generation and returns
the possible land use classes for an index of the candidate taking
into account the tabu memory.
input data:
possible_land_use is an array with possible land use classes
derived from the transition matrix
index is the index which should be checked
new_cand is a part of a new candidate
"""
# array with segments of impossible candidates (tabu memory)
global impossible_cand
copy_poss_el = np.copy(possible_land_use)
copy_new_cand = new_cand[:]
# create all new candidate parts which can be created with one of the possible_elements and the new_cand
# and check if the new one is excluded by the tabu memory
for element in possible_land_use:
copy_new_cand[index] = element
if copy_new_cand in impossible_cand:
copy_poss_el = np.delete(copy_poss_el,np.nonzero(copy_poss_el == element)[0][0])
return copy_poss_el
#------------------------------------------------------------------------------
# Check tabu memory for redundant entries
#------------------------------------------------------------------------------
def check_redundancy(new_cand):
"""Return tabu memory with new_cand as additional element
but without redundant entries.
example: if [1,1,3,4] is excluded and now [1,1,3] too, than [1,1,3,4] is
not relevant anymore
input data:
new_cand is a part of an excluded candidate
impossible_cand is the tabu memory
"""
# dictionary for possible land use options per element
global possible_elements
# array with segments of impossible candidates (tabu memory)
global impossible_cand
# original start individual derived from input data
global start_individual
search_object = new_cand[:]
sublist = []
last_elements = []
# check if new_cand is part of the impossible_cand
# if True -> do nothing
# if False -> check for redundant entries and add new_cand
if not (check_impossible_candidates(search_object)):
# check all items with one element more than the part of the new candidate
# for redundancy
if len(search_object) != len(start_individual):
for item in impossible_cand:
check_var = True
# check items with one more element than the new_cand
if len(item) == (len(search_object)+1):
# check if all elements except the last one are equal with the search_object
for index in range(len(search_object)):
if search_object[index] != item[index]:
check_var = False
break
# save potential redundant items and the last elements in separate lists
if check_var == True:
sublist.append(item)
last_elements.append(item[len(item)-1])
# delete redundant items
if cfg.mapConfig.file_transformation != 'None':