CoMOLA is a free Python tool to optimize the allocation of land use for multiple objectives. It builds upon the open source "inspyred" Python library and includes functions for reading, encoding and writing land use maps as well as genome generation and repair mutation algorithms for considering constraints during the optimization procedure. It runs on Windows and Linux and allows for the integration of any model whose prediction (e.g. a value for an ecosystem service) is based on a land use raster map. In its basic form, CoMOLA can be used immediately by inputting a raster map representing the status-quo land use, ready-to-run models written in R including their input data, and (optional) information on constraints. As constraints, the tool can consider (1) transition rules defining which type of land use can be converted into which other type and (2) minimum and maximum area proportions of each land use type within the study area. All relevant settings, such as paths to input data and models as well as optimization-specific parameters (e.g. population size, crossover and mutation rates) and settings related to constraint-handling and raster map-analysis are managed in one single control file ("config.ini").

The tool was systematically tested for different levels of complexity (Strauch et al., 2019) and applied for agricultural case studies in the Netherlands (Kromme Rijn, 219 km²) and Central Germany (Lossa, 141 km²) to best reconcile food production, biodiversity and other ecosystem services (Verhagen et al., 2018; Kaim et al., 2020). Witing et al. (2021) used CoMOLA to optimize riparian reforestation along the Zwalm River (Belgium) and Schwarz et al. (2020) optimized a virtual urban region to foster plant species richness, climate regulation and compactness. Bartkowski et al. (2020) used the tool to investigate the applicability of aligning Agent-Based Modeling with multi-objective land-use allocation. Currently, CoMOLA is providing the optimization environment for the EU research project OPTAIN (grant agreement No. 862756) to explore optimal strategies to retain water and nutrients in 14 small agricultural catchments across Europe.

CoMOLA was developed and tested for Python 2.7.

• Required Python packages
  • matplotlib
  • numpy
  • pylab

Furthermore you need to install R to run external models.

(see example files in input folder)

Land use map (required)

Provide a land use raster map in ascii format (with consecutive integer values representing the different land use classes, starting with value 1). If no user-defined patch ID map is given, CoMOLA generates a patch ID map where neighboring raster cells of the same type are aggregated.

Example (land_use.asc):

ncols         10
nrows         10
xllcorner     4376461.4080843
yllcorner     5553063.2189224
cellsize      75
NODATA_value  -2
2 6 6 1 1 2 2 3 4 4
4 6 6 1 5 5 3 6 7 8
4 2 1 1 4 4 3 3 5 5
5 2 6 6 8 8 1 1 1 1
5 7 7 2 5 8 8 4 3 3
1 7 7 4 6 8 8 4 3 1
5 2 2 2 6 8 7 5 3 1
5 5 6 6 3 4 7 5 1 1
2 2 6 6 3 4 7 3 3 3
2 6 7 7 3 2 2 1 1 1

Automatically generated patch ID map (if not pre-defined by user):

ncols         10
nrows         10
xllcorner     4376461.4080843
yllcorner     5553063.2189224
cellsize      75
NODATA_value  -2
1 2 2 3 3 4 4 5 6 6
7 2 2 3 8 8 9 10 11 0
7 12 3 3 13 13 9 9 14 14
15 12 16 16 0 0 17 17 17 17
15 18 18 19 20 0 0 21 22 22
23 18 18 24 25 0 0 21 22 26
27 28 28 28 25 0 29 30 22 26
27 27 31 31 32 33 29 30 26 26
34 34 31 31 32 33 29 35 35 35
34 36 37 37 32 38 38 39 39 39

The patch ID map is then encoded as a string of integers (each value is called a gene, representing the land use of a patch) to form the genome of the start individual, i.e. the first individual of the initial population.

Pre-defined patch ID map (optional)

If appropriate provide your own patch ID map as ascii file to delineate the spatial optimization units as needed in your specific case (with consecutive integer values representing the different patches, starting with value 1). For a cell-level optimization, an individual ID must be assigned to each cell. The spatial resolution must be the same as for the land use map. A value of 0 defines patches with a static land use. Whether a land use class is static or not must be also defined in the transition matrix (see below). Static land use classes will be excluded from the optimization.

Example (patch_IDmap_eachcell_constraint.asc):

ncols         10
nrows         10
xllcorner     4376461.4080843
yllcorner     5553063.2189224
cellsize      75
NODATA_value  -2
1 2 3 4 5 6 7 8 9 10
11 12 13 14 15 16 17 18 19 0
20 21 22 23 24 25 26 27 28 29
30 31 32 33 0 0 34 35 36 37
38 39 40 41 42 0 0 43 44 45
46 47 48 49 50 0 0 51 52 53
54 55 56 57 58 0 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

CoMOLA can handle two types of land use change constraints (simultaneously or stand-alone):

(1) Transition matrix

A matrix provided as .txt file defining which type of land use (given in rows) can be converted into which other type (given in columns). 1 = transition is possible, 0 = transition is not possible.

Example (transition_matrix.txt):

-2 1 2 3 4 5 6 7 8
1 1 1 1 1 1 1 1 0
2 1 1 1 1 1 1 1 0
3 1 1 1 1 1 1 1 0
4 1 1 1 1 1 1 1 0
5 1 1 1 1 1 1 1 0
6 0 0 0 0 0 1 1 0
7 0 0 0 0 0 0 1 0
8 0 0 0 0 0 0 0 1

(2) Total area

A table provided as .txt file defining minimum and maximum area proportions of each land use type as percentage on total area.

Example (min_max.txt):

land_use 1 2 3 4 5 6 7 8
min 0 0 0 0 0 10 10 0
max 100 100 100 100 100 25 30 100

CoMOLA handles up to four external models which must be provided as R scripts and stored in separate directories within the models folder (together with their specific input data). Each model evaluates the land use ascii map (as given above) for a specific objective that should be maximized (e.g. a certain ecosystem service). The output of each model is a single value representative for the whole study area (e.g. total agricultural yield) and needs to be written in a .csv file. If the objective value should be minimzed during optimization, multiply the model output with -1.

All relevant settings, such as paths to input data and models as well as optimization-specific parameters and settings related to constraint-handling and raster map-analysis are managed in one single control file called "config.ini". Lines starting with a semicolon (;) are commented out and not used by the algorithm.

config.ini example and description of variables:

[Dictionary]
; -----------------------------------------
; config_model
, Variable                     Description:
; ----------                   ------------
; file_path_R                | file path for R
; file_path_Python           | file path for Python
; modelx_folder              | file name of model x folder (1 <= x <= 4)
; file_modelx                | file name of the model x script
; file_outputx               | file name of the output file from model x
; update_filesx              | file names which files of the model folder x
;                              should be updated in the helping folders
; max_range                  | maximum number of land use classes
; opt_algorithm (string)     | definition of the optimization algorithm,
;                              available choices are GA or NSGA2 (default)

; RPy_available (string)     | if RPy2 is available than True or False (default)
; map                        | if True individuals are printed as ascii map files into the
;                              model folders else (default) as vectors in a csv file
; del_help_folders           | if True (default) delete and create all helping folders
;                              each time the tool starts, if False you can alternatively
;                              use the update_filex entries for updating important files
; -----------------------------------------
; config_optimization_algorithm
, Variable                     Description [default value]:
; ----------                   ------------
; pop_size                   | number of individuals per generation [100]
; max_generations            | maximum number of generations [1]
; mutation_rate              | probability for mutation [0.1]
; crossover_rate             | probability for cross over [1.0]
; priority                   | land use from NSGA2 candidate is preferred within
;                              repair mutation [True]
; maximize                   | direction of optimization [True]
; extreme_seeds              | generate extreme (but feasible) individuals for the first
;                              generation [False]
; max_repair_trials          | maximum number of repair trials within repair mutation
;                              [10000]

; terminator                 | termination criterion, see inspyred docu [default_termination]
; variator                   | variation method, see inspyred docu [default_variation]
;                              for constraint-handling and tabu-memory use repair_mutation
; selector                   | selection method, see inspyred docu [default_selection]
;                              use constrained_tournament_selection as alternative to
;                              repair_mutation
;feasible_first_pop          | if True create feasible individuals for first population
;                              [False]
;penalty_function            | 1 or 2 (only for constrained_tournament_selection)
;                                1: absolute violation measure
;                                2: normalized violation measure (default)
;plot_results                | if True plot results into a .png file [False]
; -----------------------------------------
; config_map_analysis
; Variable                     Description [default value]:
; ----------                   ------------
; file_landuse_map           | file name of land use map [none]
; four_neighbours            | analysis of four (True) or eight (False) cell neighbours
;                            | to generate patches [False]
; file_patch_map             | file name of patch ID map [none]
; file_transition            | file name of transition matrix [none]
; file_area                  | file name of total min-max area table [none]
; file_worst_fitness         | file name of worst fitness values list [none]

[config_model]
file_path_R = C:/Program Files/R/R-3.3.1/bin/R.exe
file_path_Python = C:/Python27/python.exe
model1_folder = HabStruct
file_model1 = HabStruct.R
file_output1 = HabStruct_output.csv
model2_folder = SYM
file_model2 = SYM.R
file_output2 = SYM_output.csv
model3_folder = WYLD
file_model3 = WYLD.R
file_output3 = WYLD_output.csv
;model4_folder = SAR
;file_model4 = SAR.R
;file_output4 = SAR_output.csv
max_range = 8
opt_algorithm = NSGA2

RPy2_available = False
map = True
del_help_folders = True
[config_optimization_algorithm]
pop_size = 10
max_generations = 2
mutation_rate = 0.01
crossover_rate = 0.9
priority = True
maximize = True
extreme_seeds = False
max_repair_trials = 10000
terminator = special_termination,generation_termination,diversity_termination
variator = n_point_crossover, random_reset_mutation, repair_mutation
;variator = n_point_crossover, random_reset_mutation
;selector = constrained_tournament_selection
feasible_first_pop = True
;penalty_function = 2
;plot_results = True
[config_map_analysis]
file_landuse_map = land_use.asc
four_neighbours = True
;file_patch_map = patch_IDmap_eachcell_constraint.asc
file_transition = transition_matrix.txt
file_area = min_max.txt
file_worst_fitness = worst_fitness_values_maximize.txt 

Copy the whole file structure of CoMOLA to your working directory and run __init__.py from the console:

python __init__.py 

You can limit the maximum number of threads to be used for parallel computation by adding "-t x" to the command, where x is the maximum number of threads, e.g.

python __init__.py -t 2 

Once CoMOLA has been started, a log file (Time_Date_optimization_log.txt) is generated in the output folder documenting the process of optimization.

A successful run of CoMOLA will provide the following outputs:

• fitness values for the best solutions (written at the end of the log file and in (best_solutions.csv)
• ascii maps for the best solutions
• the genome and fitness values of all individuals tested in the optimization (individuals_file.csv)

You may use the R-script (CoMOLA_postprocessing.R) provided in folder output_analysis to extract, evaluate and plot the best solutions.

Example plot: