diff --git a/.github/workflows/ci-compiled.yml b/.github/workflows/ci-compiled.yml
index 70f20bf..e8806e7 100644
--- a/.github/workflows/ci-compiled.yml
+++ b/.github/workflows/ci-compiled.yml
@@ -20,6 +20,10 @@ jobs:
         sudo apt-get install -y wget git gawk findutils
         xargs -a <(awk '! /^ *(#|$)/' ".github/workflows/apt.txt") -r -- \
             sudo apt-get install -y --no-install-recommends --no-install-suggests
+    - name: Get Python dependencies
+      run: |
+        python -m pip install --upgrade pip
+        pip install -r ./.github/workflows/requirements.txt
     - name: Set up R
       uses: r-lib/actions/setup-r@v1
       with:
diff --git a/.github/workflows/ci-ppa.yml b/.github/workflows/ci-ppa.yml
index 4b416ec..1720f97 100644
--- a/.github/workflows/ci-ppa.yml
+++ b/.github/workflows/ci-ppa.yml
@@ -23,6 +23,10 @@ jobs:
         sudo apt-get install -y -qq \
             grass grass-dev grass-doc \
             gawk findutils
+    - name: Get Python dependencies
+      run: |
+        python -m pip install --upgrade pip
+        pip install -r ./.github/workflows/requirements.txt
     - name: Set up R
       uses: r-lib/actions/setup-r@v1
       with:
diff --git a/.github/workflows/requirements.txt b/.github/workflows/requirements.txt
new file mode 100644
index 0000000..fb6c7ed
--- /dev/null
+++ b/.github/workflows/requirements.txt
@@ -0,0 +1 @@
+pandas
diff --git a/r.futures/r.futures.simulation/inputs.c b/r.futures/r.futures.simulation/inputs.c
index 3bf7606..b653950 100644
--- a/r.futures/r.futures.simulation/inputs.c
+++ b/r.futures/r.futures.simulation/inputs.c
@@ -35,7 +35,7 @@
 int get_developed_val_from_step(int step, bool abandon)
 {
     if (abandon)
-        return -step - 2;
+        return -step - 1;
     else
         return step + 1;
 }
@@ -104,10 +104,10 @@ void initialize_incentive(struct Potential *potential_info, float exponent)
  */
 void read_input_rasters(struct RasterInputs inputs, struct Segments *segments,
                         struct SegmentMemory segment_info, map_int_t *region_map,
-                        map_int_t *reverse_region_map,
+                        map_int_t *reverse_region_map, map_int_t *internal_region_map,
                         map_int_t *potential_region_map,
                         map_int_t *HUC_map, map_float_t *max_flood_probability_map,
-                        map_int_t *DDF_region_map)
+                        map_int_t *DDF_region_map, bool steering)
 {
     int row, col;
     int rows, cols;
@@ -275,10 +275,22 @@ void read_input_rasters(struct RasterInputs inputs, struct Segments *segments,
         for (col = 0; col < cols; col++) {
             isnull = false;
             /* developed */
-            /* undeveloped 0 -> -1, developed 1 -> 0 */
             if (!Rast_is_null_value(&((CELL *) developed_row)[col], CELL_TYPE)) {
                 c = ((CELL *) developed_row)[col];
-                ((CELL *) developed_row)[col] = c - 1;
+                if (steering) {
+                    if (!segments->use_climate) {
+                        /* undeveloped -1 -> -1000 */
+                        if (c == -1)
+                            ((CELL *) developed_row)[col] = DEV_TYPE_UNDEVELOPED;
+                    }
+                }
+                /* undeveloped 0 -> -1000, developed 1 -> 0 */
+                else {
+                    if (c == 0)
+                        ((CELL *) developed_row)[col] = DEV_TYPE_UNDEVELOPED;
+                    else
+                        ((CELL *) developed_row)[col] = c - 1;
+                }
             }
             else
                 isnull = true;
@@ -288,6 +300,7 @@ void read_input_rasters(struct RasterInputs inputs, struct Segments *segments,
                 region_pindex = map_get_int(region_map, c);
                 if (!region_pindex) {
                     map_set_int(region_map, c, count_regions);
+                    map_set_int(internal_region_map, c, count_regions);
                     map_set_int(reverse_region_map, count_regions, c);
                     region_index = count_regions;
                     count_regions++;
@@ -568,11 +581,12 @@ void read_predictors(struct RasterInputs inputs, struct Segments *segments,
 
 
 static int _read_demand_file(FILE *fp, const char *separator,
-                             float **table, int *demand_years,
-                             map_int_t *region_map)
+                             float ***table, struct ilist **header,
+                             int n_years, int *demand_years,
+                             map_int_t *region_map, map_int_t *reverse_region_map)
 {
 
-    size_t buflen = 4000;
+    size_t buflen = 35000;
     char buf[buflen];
     // read in the first row (header)
     if (G_getl2(buf, buflen, fp) == 0)
@@ -589,13 +603,27 @@ static int _read_demand_file(FILE *fp, const char *separator,
     if (ntokens == 0)
         G_fatal_error("No columns in the header row");
 
-    struct ilist *ids = G_new_ilist();
+    *header = G_new_ilist();
     int count;
     // skip first column which does not contain id of the region
     unsigned i;
+    int count_regions = map_nitems(region_map);
     for (i = 1; i < ntokens; i++) {
         G_chop(tokens[i]);
-        G_ilist_add(ids, atoi(tokens[i]));
+        int region_ID = atoi(tokens[i]);
+        G_ilist_add(*header, region_ID);
+        /* include subregions outside of subregion map too */
+        int *idx = map_get_int(region_map, region_ID);
+        if (!idx) {
+            map_set_int(region_map, region_ID, count_regions);
+            map_set_int(reverse_region_map, count_regions, region_ID);
+            count_regions++;
+        }
+    }
+    /* once we know how big the table is, allocate it */
+    *table = (float **) G_malloc(map_nitems(region_map) * sizeof(float *));
+    for (unsigned i = 0; i < map_nitems(region_map); i++) {
+        (*table)[i] = (float *) G_malloc(n_years * sizeof(float));
     }
 
     int years = 0;
@@ -613,22 +641,22 @@ static int _read_demand_file(FILE *fp, const char *separator,
         demand_years[years] = atoi(tokens[0]);
         for (unsigned i = 1; i < ntokens; i++) {
             // skip first column which is the year which we ignore
-            int *idx = map_get_int(region_map, ids->value[count]);
+            int *idx = map_get_int(region_map, (*header)->value[count]);
             if (idx) {
                 G_chop(tokens[i]);
-                table[*idx][years] = atof(tokens[i]);
+                (*table)[*idx][years] = atof(tokens[i]);
             }
             count++;
         }
         // each line is a year
         years++;
     }
-    G_free_ilist(ids);
     G_free_tokens(tokens);
     return years;
 }
 
-void read_demand_file(struct Demand *demandInfo, map_int_t *region_map)
+void read_demand_file(struct Demand *demandInfo, map_int_t *region_map,
+                      map_int_t *reverse_region_map)
 {
     FILE *fp_cell = NULL, *fp_population = NULL;
     if ((fp_cell = fopen(demandInfo->cells_filename, "r")) == NULL)
@@ -658,13 +686,9 @@ void read_demand_file(struct Demand *demandInfo, map_int_t *region_map)
     }
 
     demandInfo->years = (int *) G_malloc(countlines * sizeof(int));
-    demandInfo->cells_table = (float **) G_malloc(map_nitems(region_map) * sizeof(float *));
-    for (unsigned i = 0; i < map_nitems(region_map); i++) {
-        demandInfo->cells_table[i] = (float *) G_malloc(countlines * sizeof(float));
-    }
     int num_years = _read_demand_file(fp_cell, demandInfo->separator,
-                                      demandInfo->cells_table,
-                                      demandInfo->years, region_map);
+                                      &(demandInfo->cells_table), &(demandInfo->cells_header),
+                                      countlines, demandInfo->years, region_map, reverse_region_map);
     demandInfo->max_subregions = map_nitems(region_map);
     demandInfo->max_steps = num_years;
     G_verbose_message("Number of steps in area demand file: %d", num_years);
@@ -672,13 +696,9 @@ void read_demand_file(struct Demand *demandInfo, map_int_t *region_map)
 
     if (demandInfo->has_population) {
         int *years2 = (int *) G_malloc(countlines * sizeof(int));
-        demandInfo->population_table = (float **) G_malloc(map_nitems(region_map) * sizeof(float *));
-        for (unsigned i = 0; i < map_nitems(region_map); i++) {
-            demandInfo->population_table[i] = (float *) G_malloc(countlines * sizeof(float));
-        }
         int num_years2 = _read_demand_file(fp_population, demandInfo->separator,
-                                          demandInfo->population_table,
-                                          years2, region_map);
+                                           &(demandInfo->population_table), &(demandInfo->population_header),
+                                           countlines, years2, region_map, reverse_region_map);
         // check files for consistency
         if (num_years != num_years2)
             G_fatal_error(_("Area and population demand files (<%s> and <%s>) "
@@ -1164,6 +1184,8 @@ void read_flood_file(struct FloodInputs *flood_inputs)
     /* read to get length and checks */
     i = 0;
     while (G_getl2(buf, buflen, fp)) {
+        if (buf[0] == '\0')
+            continue;
         // process each column in row
         tokens = G_tokenize2(buf, flood_inputs->separator, td);
         ntokens = G_number_of_tokens(tokens);
diff --git a/r.futures/r.futures.simulation/inputs.h b/r.futures/r.futures.simulation/inputs.h
index 6f9d462..b100ec0 100644
--- a/r.futures/r.futures.simulation/inputs.h
+++ b/r.futures/r.futures.simulation/inputs.h
@@ -3,10 +3,11 @@
 
 #include <stdbool.h>
 #include <grass/segment.h>
+#include <grass/gis.h>
 
 #include "map.h"
 enum development_type {DEV_TYPE_INITIAL = 0,
-                       DEV_TYPE_UNDEVELOPED = -1};
+                       DEV_TYPE_UNDEVELOPED = -1000};
 
 enum DDF_subregions_source {DDF_DEFAULT = 0,
                             DDF_POTENTIAL = -1,
@@ -16,8 +17,12 @@ struct Demand
 {
     const char *cells_filename;
     const char *population_filename;
+    const char *cells_output_filename;
+    const char *population_output_filename;
     float **cells_table;
     float **population_table;
+    struct ilist *cells_header;
+    struct ilist *population_header;
     int *years;
     int max_subregions;
     int max_steps;
@@ -172,14 +177,15 @@ size_t estimate_undev_size(struct RasterInputs inputs);
 void initialize_incentive(struct Potential *potential_info, float exponent);
 void read_input_rasters(struct RasterInputs inputs, struct Segments *segments,
                         struct SegmentMemory segment_info, map_int_t *region_map,
-                        map_int_t *reverse_region_map,
+                        map_int_t *reverse_region_map, map_int_t *internal_region_map,
                         map_int_t *potential_region_map,
                         map_int_t *HUC_map, map_float_t *max_flood_probability_map,
-                        map_int_t *DDF_region_map);
+                        map_int_t *DDF_region_map,
+                        bool steering);
 void read_predictors(struct RasterInputs inputs, struct Segments *segments,
                      const struct Potential *potential,
                      const struct SegmentMemory segment_info);
-void read_demand_file(struct Demand *demandInfo, map_int_t *region_map);
+void read_demand_file(struct Demand *demandInfo, map_int_t *region_map, map_int_t *reverse_region_map);
 void fill_predictor_map(struct RasterInputs inputs, map_int_t *predictor_map, int num_predictors);
 void read_potential_file(struct Potential *potentialInfo, map_int_t *region_map,
                          map_int_t *predictor_map);
diff --git a/r.futures/r.futures.simulation/main.c b/r.futures/r.futures.simulation/main.c
index 8c1826f..0a89b7c 100644
--- a/r.futures/r.futures.simulation/main.c
+++ b/r.futures/r.futures.simulation/main.c
@@ -157,20 +157,26 @@ int main(int argc, char **argv)
                 *cellDemandFile, *populationDemandFile, *separator,
                 *density, *densityCapacity, *outputDensity, *redevelopmentLag,
                 *redevelopmentPotentialFile, *redistributionMatrix, *redistributionMatrixOutput,
+                *redistributionMatrixExternalOutput,
                 *HAND, *HAND_percentile, *floodInputFile, *floodLog,
                 *depthDamageFunc, *ddf_subregions, *response_func, *responseStddev,
                 *adaptations, *adaptiveCapacity, *HUCs, *outputAdaptation,
-                *patchFile, *numSteps, *output, *outputSeries, *seed, *memory;
+                *patchFile, *numSteps, *output, *outputSeries,
+                *steering_step, *outputDevPressure, *outputCellDemandFile,
+                *seed, *memory;
     } opt;
 
     struct
     {
-        struct Flag *generateSeed;
+        struct Flag *generateSeed, *steering;
     } flg;
 
     int i;
     int num_predictors;
     int num_steps;
+    int max_steps;
+    int steering_step;
+    int internal_step;
     int nseg;
     unsigned region;
     int *region_id;
@@ -183,6 +189,7 @@ int main(int argc, char **argv)
     struct RasterInputs raster_inputs;
     map_int_t region_map;
     map_int_t reverse_region_map;
+    map_int_t internal_region_map;
     map_int_t potential_region_map;
     map_int_t predictor_map;
     map_float_t max_flood_probability_map;
@@ -344,6 +351,13 @@ int main(int argc, char **argv)
         _("Basename for raster maps of development generated after each step");
     opt.outputSeries->guisection = _("Output");
 
+    opt.outputDevPressure = G_define_standard_option(G_OPT_R_OUTPUT);
+    opt.outputDevPressure->key = "output_development_pressure";
+    opt.outputDevPressure->required = NO;
+    opt.outputDevPressure->label =
+        _("Output development pressure raster");
+    opt.outputDevPressure->guisection = _("Steering");
+
     opt.outputDensity = G_define_standard_option(G_OPT_R_BASENAME_OUTPUT);
     opt.outputDensity->key = "output_density";
     opt.outputDensity->required = NO;
@@ -375,6 +389,13 @@ int main(int argc, char **argv)
             _("CSV file with population size to accommodate");
     opt.populationDemandFile->guisection = _("Demand");
 
+    opt.outputCellDemandFile = G_define_standard_option(G_OPT_F_OUTPUT);
+    opt.outputCellDemandFile->key = "output_demand";
+    opt.outputCellDemandFile->required = NO;
+    opt.outputCellDemandFile->description =
+            _("Output CSV file with number of cells to convert for each step and subregion");
+    opt.outputCellDemandFile->guisection = _("Steering");
+
     opt.separator = G_define_standard_option(G_OPT_F_SEP);
     opt.separator->answer = "comma";
     opt.separator->description =
@@ -399,6 +420,12 @@ int main(int argc, char **argv)
     opt.redistributionMatrixOutput->description = _("Base name for output file containing matrix of pixels moved from one subregion to another");
     opt.redistributionMatrixOutput->guisection = _("Climate scenarios");
 
+    opt.redistributionMatrixExternalOutput = G_define_standard_option(G_OPT_F_OUTPUT);
+    opt.redistributionMatrixExternalOutput->key = "redistribution_external_output";
+    opt.redistributionMatrixExternalOutput->required = NO;
+    opt.redistributionMatrixExternalOutput->description = _("Base name for output file containing matrix of pixels moved from one subregion to another");
+    opt.redistributionMatrixExternalOutput->guisection = _("Steering");
+
     opt.HAND = G_define_standard_option(G_OPT_R_INPUT);
     opt.HAND->key = "hand";
     opt.HAND->required = NO;
@@ -537,6 +564,15 @@ int main(int argc, char **argv)
         _("Number of steps to be simulated");
     opt.numSteps->guisection = _("Basic input");
 
+    opt.steering_step = G_define_option();
+    opt.steering_step->key = "steering_step";
+    opt.steering_step->type = TYPE_INTEGER;
+    opt.steering_step->required = NO;
+    opt.steering_step->options = "0-1000";
+    opt.steering_step->description =
+        _("Steering step");
+    opt.steering_step->guisection = _("Steering");
+
     opt.potentialWeight = G_define_standard_option(G_OPT_R_INPUT);
     opt.potentialWeight->key = "potential_weight";
     opt.potentialWeight->required = NO;
@@ -584,6 +620,11 @@ int main(int argc, char **argv)
     opt.memory->required = NO;
     opt.memory->description = _("Memory in GB");
 
+    flg.steering = G_define_flag();
+    flg.steering->key = 'r';
+    flg.steering->description = _("Use steering");
+    flg.steering->guisection = _("Steering");
+
     // TODO: add mutually exclusive?
     // TODO: add flags or options to control values in series and final rasters
 
@@ -599,6 +640,8 @@ int main(int argc, char **argv)
     G_option_requires(opt.outputAdaptation, opt.adaptiveCapacity, NULL);
     G_option_requires(opt.HAND, opt.HAND_percentile, NULL);
     G_option_requires(opt.floodLog, opt.floodInputFile, NULL);
+    G_option_requires_all(flg.steering, opt.outputDevPressure,
+                          opt.outputCellDemandFile, opt.steering_step, NULL);
     if (G_parser(argc, argv))
         exit(EXIT_FAILURE);
 
@@ -642,11 +685,12 @@ int main(int argc, char **argv)
     patch_info.strategy = SKIP;
     if (opt.redevelopmentLag->answer)
         patch_info.redevelopment_lag = atoi(opt.redevelopmentLag->answer);
-    
-    num_steps = 0;
-    if (opt.numSteps->answer)
-        num_steps = atoi(opt.numSteps->answer);
-    
+
+    step = 0;
+    if (opt.steering_step->answer) {
+        steering_step = atoi(opt.steering_step->answer);
+        step = steering_step;
+    }
     num_predictors = 0;
     for (i = 0; opt.predictors->answers[i]; i++)
         num_predictors++;
@@ -765,14 +809,15 @@ int main(int argc, char **argv)
     //    read Subregions layer
     map_init(&region_map);
     map_init(&reverse_region_map);
+    map_init(&internal_region_map);
     map_init(&potential_region_map);
     map_init(&predictor_map);
     map_init(&DDF_region_map);
     G_verbose_message("Reading input rasters...");
     read_input_rasters(raster_inputs, &segments, segment_info, &region_map,
-                       &reverse_region_map, &potential_region_map,
+                       &reverse_region_map, &internal_region_map, &potential_region_map,
                        &HUC_map, &max_flood_probability_map,
-                       &DDF_region_map);
+                       &DDF_region_map, flg.steering->answer);
     if (flood_inputs.array) {
         create_bboxes(&segments.HUC, &segments.developed, &bboxes);
     }
@@ -812,15 +857,29 @@ int main(int argc, char **argv)
         demand_info.has_population = true;
         demand_info.population_filename = opt.populationDemandFile->answer;
     }
+    if (opt.outputCellDemandFile->answer)
+        demand_info.cells_output_filename = opt.outputCellDemandFile->answer;
     demand_info.separator = G_option_to_separator(opt.separator);
-    read_demand_file(&demand_info, &region_map);
-    if (num_steps == 0)
-        num_steps = demand_info.max_steps;
+    read_demand_file(&demand_info, &region_map, &reverse_region_map);
+
+    max_steps = demand_info.max_steps;
+    internal_step = 0;
+    num_steps = max_steps;
+    if (opt.numSteps->answer)
+        num_steps = atoi(opt.numSteps->answer);
+
 
     /* check redistribution matrix output files */
     if (opt.redistributionMatrixOutput->answer) {
         redistr_matrix.output_basename = opt.redistributionMatrixOutput->answer;
-        if (check_matrix_filenames_exist(&redistr_matrix, num_steps)) {
+        if (check_matrix_filenames_exist(&redistr_matrix, false, max_steps)) {
+            if (!G_check_overwrite(argc, argv))
+                G_fatal_error(_("At least one of the requested matrix output files exists. Use --o to overwrite."));
+        }
+    }
+    if (opt.redistributionMatrixExternalOutput->answer) {
+        redistr_matrix.output_external_basename = opt.redistributionMatrixExternalOutput->answer;
+        if (check_matrix_filenames_exist(&redistr_matrix, true, max_steps)) {
             if (!G_check_overwrite(argc, argv))
                 G_fatal_error(_("At least one of the requested matrix output files exists. Use --o to overwrite."));
         }
@@ -848,14 +907,15 @@ int main(int argc, char **argv)
     /* read Patch sizes file */
     G_verbose_message("Reading patch size file...");
     patch_sizes.filename = opt.patchFile->answer;
-    read_patch_sizes(&patch_sizes, &region_map, discount_factor);
+    read_patch_sizes(&patch_sizes, &internal_region_map, discount_factor);
 
-    undev_cells = initialize_developables(map_nitems(&region_map), undev_estimate);
+    undev_cells = initialize_developables(map_nitems(&internal_region_map), undev_estimate);
+    /* all demand related arrays need to account for external subregions, hence region_map */
     patch_overflow = G_calloc(map_nitems(&region_map), sizeof(int));
 
     /* redevelopment */
     if (segments.use_density) {
-        dev_cells = initialize_developables(map_nitems(&region_map), undev_estimate);
+        dev_cells = initialize_developables(map_nitems(&internal_region_map), undev_estimate);
         population_overflow = G_calloc(map_nitems(&region_map), sizeof(float));
     }
     else {
@@ -866,17 +926,18 @@ int main(int argc, char **argv)
     overgrow = true;
     G_verbose_message("Starting simulation...");
     leaving_population = 0;
-    for (step = 0; step < num_steps; step++) {
+    for (; step < max_steps && internal_step < num_steps; step++, internal_step++) {
         recompute_probabilities(undev_cells, &segments, &potential_info, false);
         if (segments.use_density)
             recompute_probabilities(dev_cells, &segments, &redev_potential_info, true);
-        if (step == num_steps - 1)
+        if (step == max_steps - 1)
             overgrow = false;
-        for (region = 0; region < map_nitems(&region_map); region++) {
+        for (region = 0; region < map_nitems(&internal_region_map); region++) {
+            /* Indices of regions present in subregions always start with 0 because they are read in first */
             region_id = map_get_int(&reverse_region_map, region);
             G_verbose_message("Computing step %d (out of %d), region %d (%d out of %d)",
-                              step + 1, num_steps, *region_id,
-                              region + 1, map_nitems(&region_map));
+                              step + 1, max_steps, *region_id,
+                              region + 1, map_nitems(&internal_region_map));
             compute_step(undev_cells, dev_cells, &demand_info, search_alg, &segments,
                          &patch_sizes, &patch_info, &devpressure_info, patch_overflow,
                          population_overflow, &redistr_matrix, step, region, &reverse_region_map, overgrow);
@@ -889,37 +950,45 @@ int main(int argc, char **argv)
                 update_flood_probability(step, &flood_inputs, &segments, &HUC_map, &max_flood_probability_map);
             for (HUC = 0; HUC < map_nitems(&HUC_map); HUC++)
                 climate_step(&segments, &demand_info, &bboxes,
-                             &redistr_matrix, &region_map, &reverse_region_map,
+                             &redistr_matrix, &region_map, &reverse_region_map, &internal_region_map,
                              step, &leaving_population,
                              &HUC_bbox_vals, HAND_percentile,
                              &max_flood_probability_map, &flood_inputs, &flood_log, &DDF,
                              &response_relation, HUC);
             if (opt.redistributionMatrixOutput->answer)
-                write_redistribution_matrix(&redistr_matrix, step, num_steps);
+                write_redistribution_matrix(&redistr_matrix, false, step, max_steps);
+            if (opt.redistributionMatrixExternalOutput->answer)
+                write_redistribution_matrix(&redistr_matrix, true, step, max_steps);
         }
         /* export developed for that step */
         if (opt.outputSeries->answer) {
-            name_step = name_for_step(opt.outputSeries->answer, step, num_steps);
+            name_step = name_for_step(opt.outputSeries->answer, step, max_steps);
             output_developed_step(&segments.developed, name_step,
-                                  demand_info.years[step], -1, num_steps,
+                                  demand_info.years[step], -1, max_steps,
                                   segments.use_climate ? true : false);
         }
         /* export density for that step */
         if (opt.outputDensity->answer) {
-            name_step = name_for_step(opt.outputDensity->answer, step, num_steps);
+            name_step = name_for_step(opt.outputDensity->answer, step, max_steps);
             output_step(&segments.density, &segments.developed, name_step, FCELL_TYPE);
         }
         /* export density for that step */
         if (opt.outputAdaptation->answer) {
-            name_step = name_for_step(opt.outputAdaptation->answer, step, num_steps);
+            name_step = name_for_step(opt.outputAdaptation->answer, step, max_steps);
             output_step(&segments.adaptation, &segments.developed, name_step, CELL_TYPE);
         }
+        if (opt.outputCellDemandFile->answer)
+            output_demand_file(&demand_info, &region_map, patch_overflow, step);
     }
 
     /* write */
     output_developed_step(&segments.developed, opt.output->answer,
                           demand_info.years[0], demand_info.years[step-1],
-                          num_steps, segments.use_climate ? true : false);
+                          max_steps, segments.use_climate ? true : false);
+
+    if (opt.outputDevPressure->answer)
+        output_step(&segments.devpressure, &segments.developed, opt.outputDevPressure->answer, FCELL_TYPE);
+
     if (opt.floodLog->answer)
         write_flood_log(&flood_log, opt.floodLog->answer, &HUC_map);
 
@@ -960,6 +1029,7 @@ int main(int argc, char **argv)
     }
     map_deinit(&region_map);
     map_deinit(&reverse_region_map);
+    map_deinit(&internal_region_map);
     map_deinit(&predictor_map);
     map_deinit(&potential_region_map);
     map_deinit(&DDF_region_map);
@@ -968,11 +1038,13 @@ int main(int argc, char **argv)
             G_free(demand_info.cells_table[i]);
         G_free(demand_info.cells_table);
         G_free(demand_info.years);
+        G_free_ilist(demand_info.cells_header);
     }
     if (demand_info.has_population) {
         for (int i = 0; i < demand_info.max_subregions; i++)
             G_free(demand_info.population_table[i]);
         G_free(demand_info.population_table);
+        G_free_ilist(demand_info.population_header);
     }
     if (potential_info.predictors) {
         for (int i = 0; i < potential_info.max_predictors; i++)
diff --git a/r.futures/r.futures.simulation/output.c b/r.futures/r.futures.simulation/output.c
index 47b4dd3..d50758e 100644
--- a/r.futures/r.futures.simulation/output.c
+++ b/r.futures/r.futures.simulation/output.c
@@ -23,6 +23,7 @@
 
 #include "output.h"
 #include "inputs.h"
+#include "map.h"
 
 
 static void create_timestamp(int year, struct TimeStamp* timestamp)
@@ -73,9 +74,9 @@ char *name_for_step(const char *basename, const int step, const int nsteps)
  * \param year_from year to put as timestamp
  * \param year_to if > 0 it is end year of timestamp interval
  * \param nsteps total number of steps (needed for color table)
- * \param contains_abandoned If true, represent undeveloped areas as NULLS, and
- *        shift values <= -2 (abandoned) by +1 to represent in absolute value
- *        step when they were abandoned.
+ * \param contains_abandoned If false, undeveloped areas represented internally as -1000
+ *        are changed on output to -1, otherwise -1000 is kept on output
+ *        with -1 to -nsteps representing the step in abs value when px was abandoned
  */
 void output_developed_step(SEGMENT *developed_segment, const char *name,
                            int year_from, int year_to, int nsteps, bool contains_abandoned)
@@ -103,13 +104,9 @@ void output_developed_step(SEGMENT *developed_segment, const char *name,
             if (Rast_is_c_null_value(&developed)) {
                 continue;
             }
-            if (contains_abandoned) {
-                if (developed < DEV_TYPE_UNDEVELOPED) {
-                    developed += 1;
-                }
-                else if (developed == DEV_TYPE_UNDEVELOPED) {
-                    continue;
-                }
+            if (!contains_abandoned) {
+                if (developed == DEV_TYPE_UNDEVELOPED)
+                    developed = -1;
             }
             out_row[col] = developed;
         }
@@ -122,14 +119,18 @@ void output_developed_step(SEGMENT *developed_segment, const char *name,
     // TODO: the map max is 36 for 36 steps, it is correct?
 
     if (contains_abandoned) {
+        val1 = DEV_TYPE_UNDEVELOPED;
+        val2 = DEV_TYPE_UNDEVELOPED;
+        Rast_add_c_color_rule(&val1, 180, 255, 160, &val2, 180, 255, 160,
+                              &colors);
         val1 = -nsteps;
         val2 = -1;
         Rast_add_c_color_rule(&val1, 230, 163, 255, &val2, 60, 0, 80,
                               &colors);
-     }
+    }
     else {
-        val1 = DEV_TYPE_UNDEVELOPED;
-        val2 = DEV_TYPE_UNDEVELOPED;
+        val1 = -1;
+        val2 = -1;
         Rast_add_c_color_rule(&val1, 180, 255, 160, &val2, 180, 255, 160,
                               &colors);
     }
@@ -207,3 +208,44 @@ void output_step(SEGMENT *output_segment, SEGMENT *developed_segment,
     G_message(_("Raster map <%s> created"), name);
 
 }
+
+void output_demand_file(struct Demand *demandInfo, map_int_t *region_map,
+                        int *patch_overflow, int step)
+{
+    FILE *fp_cell = NULL;
+    if ((fp_cell = fopen(demandInfo->cells_output_filename, "w")) == NULL)
+        G_fatal_error(_("Cannot open area demand file <%s> for writing"),
+                      demandInfo->cells_output_filename);
+    fprintf(fp_cell, "year");
+    /* write region header */
+    for (int region = 0; region < demandInfo->cells_header->n_values; region++) {
+        fprintf(fp_cell, ",%d", demandInfo->cells_header->value[region]);
+    }
+    fprintf(fp_cell, "\n");
+
+    int *patch_overflow_copy = G_calloc(map_nitems(region_map), sizeof(int));
+    memcpy(patch_overflow_copy, patch_overflow, sizeof(int) * map_nitems(region_map));
+    for (int row = 0; row < demandInfo->max_steps; row++) {
+        /* write year */
+        fprintf(fp_cell, "%d", demandInfo->years[row]);
+        for (int region = 0; region < demandInfo->cells_header->n_values; region++) {
+            int *region_idx = map_get_int(region_map, demandInfo->cells_header->value[region]);
+            /* resolve patch overflow */
+            int n_to_convert = lroundf(demandInfo->cells_table[*region_idx][row]);
+
+            if ((row == step + 1) && patch_overflow_copy[*region_idx] > 0) {
+                if (n_to_convert - patch_overflow_copy[*region_idx] > 0) {
+                    n_to_convert -= patch_overflow_copy[*region_idx];
+                }
+                else {
+                    patch_overflow_copy[*region_idx] -= n_to_convert;
+                    n_to_convert = 0;
+                }
+            }
+            fprintf(fp_cell, ",%d", n_to_convert);
+        }
+        fprintf(fp_cell, "\n");
+    }
+    G_free(patch_overflow_copy);
+    fclose(fp_cell);
+}
diff --git a/r.futures/r.futures.simulation/output.h b/r.futures/r.futures.simulation/output.h
index 23854c4..78f67c5 100644
--- a/r.futures/r.futures.simulation/output.h
+++ b/r.futures/r.futures.simulation/output.h
@@ -5,6 +5,9 @@
 #include <grass/segment.h>
 #include <stdbool.h>
 
+#include "map.h"
+#include "inputs.h"
+
 
 
 char *name_for_step(const char *basename, const int step, const int nsteps);
@@ -12,4 +15,5 @@ void output_developed_step(SEGMENT *developed_segment, const char *name, int yea
                            int nsteps, bool contains_abandoned);
 void output_step(SEGMENT *output_segment, SEGMENT *developed_segment,
                  const char *name, RASTER_MAP_TYPE data_type);
+void output_demand_file(struct Demand *demandInfo, map_int_t *region_map, int *patch_overflow, int step);
 #endif // FUTURES_OUTPUT_H
diff --git a/r.futures/r.futures.simulation/redistribute.c b/r.futures/r.futures.simulation/redistribute.c
index 21e0c50..953ea08 100644
--- a/r.futures/r.futures.simulation/redistribute.c
+++ b/r.futures/r.futures.simulation/redistribute.c
@@ -1,4 +1,4 @@
-/*!
+/*!
    \file redistribute.c
 
    \brief Functions to redistribute population
@@ -63,7 +63,7 @@ int pick_region(const float *probabilities, int size)
  * \param leaving_population population leaving outside of study area
  */
 void redistribute(struct RedistributionMatrix *matrix, struct Demand *demand,
-                  int regionID, int num_px, map_int_t *region_map,
+                  int regionID, int num_px, map_int_t *region_map, map_int_t *internal_region_map,
                   int step, float *leaving_population)
 {
     int to_idx;
@@ -71,6 +71,7 @@ void redistribute(struct RedistributionMatrix *matrix, struct Demand *demand,
     int *to_ID;
     int *demand_to_idx;
     int *demand_from_idx;
+    int *internal_to_idx;
     float density_from, density_to;
     float to_px;
 
@@ -97,7 +98,13 @@ void redistribute(struct RedistributionMatrix *matrix, struct Demand *demand,
             to_px = num_px * density_from / density_to;
         /* increase number of px to grow next step */
         if (step + 1 < demand->max_steps) {
-            demand->cells_table[*demand_to_idx][step + 1] += to_px;
+            internal_to_idx = map_get_int(internal_region_map, *to_ID);
+            /* if destination is inside study area, apply to demand
+               otherwise record in external matrix */
+            if (internal_to_idx)
+                demand->cells_table[*demand_to_idx][step + 1] += to_px;
+            else
+                matrix->external_px[*from_idx][to_idx] += to_px;
             G_debug(2, "%f cells moved from %d to %d in step %d", to_px, regionID, *to_ID, step);
             matrix->moved_px[*from_idx][to_idx] += to_px;
         }
@@ -162,6 +169,7 @@ void read_redistribution_matrix(struct RedistributionMatrix *matrix)
     matrix->dim_from = matrix->dim_to;
     matrix->probabilities = (float **) G_malloc(matrix->dim_from * sizeof(float *));
     matrix->moved_px = (float **) G_malloc(matrix->dim_from * sizeof(float *));
+    matrix->external_px = (float **) G_malloc(matrix->dim_from * sizeof(float *));
     matrix->max_dim_from = matrix->dim_from;
 
     row = 0;
@@ -175,6 +183,7 @@ void read_redistribution_matrix(struct RedistributionMatrix *matrix)
         /* allocate row */
         matrix->probabilities[row] = (float *) G_malloc(matrix->dim_to * sizeof(float));
         matrix->moved_px[row] = (float *) G_malloc(matrix->dim_to * sizeof(float));
+        matrix->external_px[row] = (float *) G_malloc(matrix->dim_to * sizeof(float));
         for (col = 0; col < ntokens2; col++) {
             if (col == 0) {
                 map_set(&matrix->from_map, tokens2[col], row);
@@ -184,6 +193,7 @@ void read_redistribution_matrix(struct RedistributionMatrix *matrix)
                 /* convert from % */
                 matrix->probabilities[row][col - 1] = atof(tokens2[col]) / 100.;
                 matrix->moved_px[row][col - 1] = 0;
+                matrix->external_px[row][col - 1] = 0;
             }
         }
         row++;
@@ -191,6 +201,7 @@ void read_redistribution_matrix(struct RedistributionMatrix *matrix)
             new_size = 2 * matrix->max_dim_from;
             matrix->probabilities = (float **) G_realloc(matrix->probabilities, new_size * sizeof(float *));
             matrix->moved_px = (float **) G_realloc(matrix->moved_px, new_size * sizeof(float *));
+            matrix->external_px = (float **) G_realloc(matrix->external_px, new_size * sizeof(float *));
             matrix->max_dim_from = new_size;
         }
 
@@ -210,14 +221,14 @@ void read_redistribution_matrix(struct RedistributionMatrix *matrix)
  * \param nsteps Number of steps in simulation
  * \return
  */
-static char* get_matrix_filename(const struct RedistributionMatrix *matrix, int step, int nsteps)
+static char* get_matrix_filename(const struct RedistributionMatrix *matrix, bool external, int step, int nsteps)
 {
     const char *name;
     const char *ext;
     char *filename;
 
     ext = ".csv";
-    name = name_for_step(matrix->output_basename, step, nsteps);
+    name = name_for_step(external ? matrix->output_external_basename : matrix->output_basename, step, nsteps);
     filename = G_malloc((strlen(name) + strlen(ext) + 1) * sizeof(char));
     sprintf(filename, "%s%s", name, ext);
     return filename;
@@ -229,12 +240,12 @@ static char* get_matrix_filename(const struct RedistributionMatrix *matrix, int
  * \param nsteps Number of steps in simulation
  * \return True if any of the files exist otherwise false
  */
-bool check_matrix_filenames_exist(const struct RedistributionMatrix *matrix, int nsteps)
+bool check_matrix_filenames_exist(const struct RedistributionMatrix *matrix, bool external, int nsteps)
 {
     char *filename;
     int step;
     for (step = 0; step < nsteps; step++) {
-        filename = get_matrix_filename(matrix, step, nsteps);
+        filename = get_matrix_filename(matrix, external, step, nsteps);
         if (access(filename, F_OK) != -1) {
             G_free(filename);
             return true;
@@ -247,11 +258,12 @@ bool check_matrix_filenames_exist(const struct RedistributionMatrix *matrix, int
 /*!
  * \brief Write redistribution matrix (number of moved pixels among counties)
  * \param matrix structure
+ * \param external if write external moved px (for steered distributed runs)
  * \param step step number
  * \param nsteps total number of steps (for file name padding)
  */
 void write_redistribution_matrix(struct RedistributionMatrix *matrix,
-                                 int step, int nsteps)
+                                 bool external, int step, int nsteps)
 {
     FILE *fp;
     const char *name;
@@ -261,7 +273,10 @@ void write_redistribution_matrix(struct RedistributionMatrix *matrix,
     int *ID;
 
     ext = ".csv";
-    name = name_for_step(matrix->output_basename, step, nsteps);
+    if (external)
+        name = name_for_step(matrix->output_external_basename, step, nsteps);
+    else
+        name = name_for_step(matrix->output_basename, step, nsteps);
     filename = G_malloc((strlen(name) + strlen(ext) + 1) * sizeof(char));
     sprintf(filename, "%s%s", name, ext);
     fp = fopen(filename, "w+");
@@ -271,14 +286,15 @@ void write_redistribution_matrix(struct RedistributionMatrix *matrix,
         fprintf(fp, ",%d", *ID);
     }
     fprintf(fp, "\n");
+    float **matrix_px = external ? matrix->external_px : matrix->moved_px;
     for (row = 0; row < matrix->dim_from; row++) {
         ID = map_get_int(&matrix->reverse_from_map, row);
         fprintf(fp, "%d", *ID);
         for (col = 0; col < matrix->dim_to; col++) {
-            if (matrix->moved_px[row][col] > 0)
-                fprintf(fp, ",%f", matrix->moved_px[row][col]);
+            if (matrix_px[row][col] > 0)
+                fprintf(fp, ",%f", matrix_px[row][col]);
             else
-                fprintf(fp, ",%d", (int)matrix->moved_px[row][col]);
+                fprintf(fp, ",%d", (int)matrix_px[row][col]);
         }
         fprintf(fp, "\n");
     }
@@ -292,9 +308,11 @@ void free_redistribution_matrix(struct RedistributionMatrix *matrix)
     for (i = 0; i < matrix->dim_from; i++) {
         G_free(matrix->probabilities[i]);
         G_free(matrix->moved_px[i]);
+        G_free(matrix->external_px[i]);
     }
     G_free(matrix->probabilities);
     G_free(matrix->moved_px);
+    G_free(matrix->external_px);
     map_deinit(&matrix->from_map);
     map_deinit(&matrix->to_map);
     map_deinit(&matrix->reverse_from_map);
diff --git a/r.futures/r.futures.simulation/redistribute.h b/r.futures/r.futures.simulation/redistribute.h
index 8cdac5f..8447899 100644
--- a/r.futures/r.futures.simulation/redistribute.h
+++ b/r.futures/r.futures.simulation/redistribute.h
@@ -8,8 +8,10 @@ struct RedistributionMatrix
 {
     const char *filename;
     const char *output_basename;
+    const char *output_external_basename;
     float **probabilities;
     float **moved_px;
+    float **external_px;
     int dim_from;
     int dim_to;
     int max_dim_from;
@@ -19,10 +21,10 @@ struct RedistributionMatrix
 };
 
 void redistribute(struct RedistributionMatrix *matrix, struct Demand *demand,
-                  int regionID, int num_px, map_int_t *region_map,
+                  int regionID, int num_px, map_int_t *region_map, map_int_t *internal_region_map,
                   int step, float *leaving_population);
 void read_redistribution_matrix(struct RedistributionMatrix *matrix);
-bool check_matrix_filenames_exist(const struct RedistributionMatrix *matrix, int nsteps);
-void write_redistribution_matrix(struct RedistributionMatrix *matrix, int step, int nsteps);
+bool check_matrix_filenames_exist(const struct RedistributionMatrix *matrix, bool external, int nsteps);
+void write_redistribution_matrix(struct RedistributionMatrix *matrix, bool external, int step, int nsteps);
 void free_redistribution_matrix(struct RedistributionMatrix *matrix);
 #endif // FUTURES_REDISTRIBUTE_H
diff --git a/r.futures/r.futures.simulation/simulation.c b/r.futures/r.futures.simulation/simulation.c
index f2911c5..203094f 100644
--- a/r.futures/r.futures.simulation/simulation.c
+++ b/r.futures/r.futures.simulation/simulation.c
@@ -444,6 +444,7 @@ void compute_step(struct Developables *undev_cells, struct Developables *dev_cel
 void climate_step(struct Segments *segments, struct Demand *demand,
                   struct BBoxes *bboxes, struct RedistributionMatrix *matrix,
                   map_int_t *region_map, map_int_t *reverse_region_map,
+                  map_int_t *internal_region_map,
                   int step, float *leaving_population,
                   struct HAND_bbox_values *HAND_bbox_vals,
                   float percentile,
@@ -519,7 +520,7 @@ void climate_step(struct Segments *segments, struct Demand *demand,
                             Segment_get(&segments->subregions, (void *)&region_from_idx, row, col);
                             region_from_ID = map_get_int(reverse_region_map, region_from_idx);
                             redistribute(matrix, demand, *region_from_ID, 1,
-                                         region_map, step, leaving_population);
+                                         region_map, internal_region_map, step, leaving_population);
                         }
                         else if (response == Adapt) {
                             adapt(&segments->adaptation, flood_probability, row, col);
diff --git a/r.futures/r.futures.simulation/simulation.h b/r.futures/r.futures.simulation/simulation.h
index 86b5311..5e69f3f 100644
--- a/r.futures/r.futures.simulation/simulation.h
+++ b/r.futures/r.futures.simulation/simulation.h
@@ -45,7 +45,7 @@ void compute_step(struct Developables *undev_cells, struct Developables *dev_cel
                   bool overgrow);
 void climate_step(struct Segments *segments, struct Demand *demand,
                   struct BBoxes *bboxes, struct RedistributionMatrix *matrix,
-                  map_int_t *region_map, map_int_t *reverse_region_map,
+                  map_int_t *region_map, map_int_t *reverse_region_map, map_int_t *internal_region_map,
                   int step, float *leaving_population,
                   struct HAND_bbox_values *HAND_bbox_vals, float percentile,
                   map_float_t *flood_probability_map, const struct FloodInputs *flood_inputs,
diff --git a/r.futures/r.futures.simulation/testsuite/data/flood_depth_input_steered.csv b/r.futures/r.futures.simulation/testsuite/data/flood_depth_input_steered.csv
new file mode 100644
index 0000000..37cd37c
--- /dev/null
+++ b/r.futures/r.futures.simulation/testsuite/data/flood_depth_input_steered.csv
@@ -0,0 +1,22 @@
+step,return,map
+1,0.05,flood_1
+1,0.1,flood_2
+1,0.2,flood_3
+2,0.05,flood_1
+2,0.1,flood_2
+2,0.2,flood_3
+3,0.05,flood_1
+3,0.1,flood_2
+3,0.2,flood_3
+4,0.05,flood_1
+4,0.1,flood_2
+4,0.2,flood_3
+5,0.05,flood_1
+5,0.1,flood_2
+5,0.2,flood_3
+6,0.05,flood_1
+6,0.1,flood_2
+6,0.2,flood_3
+7,0.05,flood_1
+7,0.1,flood_2
+7,0.2,flood_3
diff --git a/r.futures/r.futures.simulation/testsuite/data/result_steering.pack b/r.futures/r.futures.simulation/testsuite/data/result_steering.pack
new file mode 100644
index 0000000..22d25ab
Binary files /dev/null and b/r.futures/r.futures.simulation/testsuite/data/result_steering.pack differ
diff --git a/r.futures/r.futures.simulation/testsuite/data/result_steering_distributed.pack b/r.futures/r.futures.simulation/testsuite/data/result_steering_distributed.pack
new file mode 100644
index 0000000..73841f7
Binary files /dev/null and b/r.futures/r.futures.simulation/testsuite/data/result_steering_distributed.pack differ
diff --git a/r.futures/r.futures.simulation/testsuite/test_r_futures_simulation.py b/r.futures/r.futures.simulation/testsuite/test_r_futures_simulation.py
index 1dbf0b6..0e6c3af 100755
--- a/r.futures/r.futures.simulation/testsuite/test_r_futures_simulation.py
+++ b/r.futures/r.futures.simulation/testsuite/test_r_futures_simulation.py
@@ -10,6 +10,8 @@ class TestPGA(TestCase):
     output = "pga_output"
     result = "result"
     potential_test = "data/potential_test.csv"
+    output_demand = "data/output_demand_file.csv"
+    output_development_pressure = "output_development_pressure"
 
     @classmethod
     def setUpClass(cls):
@@ -109,6 +111,9 @@ def setUpClass(cls):
         cls.runModule(
             "r.mapcalc", expression="acapacity = rand(-100, 100) * 0.01", seed=1
         )
+        cls.runModule(
+            "r.mapcalc", expression="urban_2002_steered = if (urban_2002 == 0, -1, 0)"
+        )
 
     @classmethod
     def tearDownClass(cls):
@@ -138,14 +143,21 @@ def tearDownClass(cls):
                 "flood_probability",
                 "acapacity",
                 "basin",
+                "urban_2002_steered",
                 cls.result,
             ],
         )
         cls.del_temp_region()
         gs.try_remove(cls.potential_test)
+        gs.try_remove(cls.output_demand)
 
     def tearDown(self):
-        self.runModule("g.remove", flags="f", type="raster", name=self.output)
+        self.runModule(
+            "g.remove",
+            flags="f",
+            type="raster",
+            name=[self.output, self.output_development_pressure],
+        )
 
     def test_pga_run(self):
         """Test if results is in expected limits"""
@@ -340,6 +352,40 @@ def test_pga_flooding_DDF(self):
             output=self.output,
         )
 
+    def test_pga_run_steering(self):
+        """Test if results is in expected limits"""
+        for i in range(0, 7):
+            self.assertModule(
+                "r.futures.simulation",
+                developed="urban_2002_steered" if i == 0 else self.output,
+                development_pressure="devpressure"
+                if i == 0
+                else self.output_development_pressure,
+                compactness_mean=0.4,
+                compactness_range=0.05,
+                discount_factor=0.1,
+                patch_sizes="data/patches.txt",
+                predictors=["slope", "lakes_dist_km", "streets_dist_km"],
+                n_dev_neighbourhood=15,
+                devpot_params="data/potential.csv",
+                random_seed=1,
+                num_neighbors=4,
+                seed_search="random",
+                development_pressure_approach="gravity",
+                gamma=1.5,
+                scaling_factor=1,
+                subregions="zipcodes",
+                demand="data/demand.csv" if i == 0 else self.output_demand,
+                output=self.output,
+                flags="r",
+                output_development_pressure=self.output_development_pressure,
+                output_demand=self.output_demand,
+                num_steps=1,
+                steering_step=i,
+                verbose=True,
+                overwrite=True,
+            )
+
 
 if __name__ == "__main__":
     test()
diff --git a/r.futures/r.futures.simulation/testsuite/test_r_futures_simulation_steering.py b/r.futures/r.futures.simulation/testsuite/test_r_futures_simulation_steering.py
new file mode 100644
index 0000000..439cf00
--- /dev/null
+++ b/r.futures/r.futures.simulation/testsuite/test_r_futures_simulation_steering.py
@@ -0,0 +1,497 @@
+#!/usr/bin/env python3
+import pandas as pd
+
+from grass.gunittest.case import TestCase
+from grass.gunittest.main import test
+import grass.script as gs
+
+
+def compare_base_stats(reference, output_map):
+    def stats(input_map):
+        data = gs.read_command(
+            "r.stats", flags="cn", input=["zipcodes", input_map], separator="comma"
+        )
+        cell_list = []
+        for line in data.strip().splitlines():
+            zipcode, step, cells = line.split(",")
+            if step == "-1":
+                cell_list.append(int(cells))
+        return cell_list
+
+    reference_stats = stats(reference)
+    output_stats = stats(output_map)
+    return [
+        100 * abs(r - o) / (r if r > o else o)
+        for r, o in zip(reference_stats, output_stats)
+    ]
+
+
+class TestSteering(TestCase):
+
+    output = "pga_output"
+    output_1 = "pga_output_1"
+    output_2 = "pga_output_2"
+    result = "result"
+    result_steering = "result_steering"
+    result_steering_distributed = "result_steering_distributed"
+    output_demand = "output_demand_file.csv"
+    output_demand_1 = "output_demand_file_1.csv"
+    output_demand_2 = "output_demand_file_2.csv"
+    output_development_pressure = "output_development_pressure"
+    output_development_pressure_1 = "output_development_pressure_1"
+    output_development_pressure_2 = "output_development_pressure_2"
+    redistribution_1 = "redistribution_1"
+    redistribution_2 = "redistribution_2"
+
+    @classmethod
+    def setUpClass(cls):
+        cls.use_temp_region()
+        cls.runModule("g.region", raster="lsat7_2002_30@PERMANENT")
+        cls.runModule("r.unpack", input="data/result.pack", output=cls.result)
+        cls.runModule(
+            "r.unpack", input="data/result_steering.pack", output=cls.result_steering
+        )
+        cls.runModule(
+            "r.unpack",
+            input="data/result_steering_distributed.pack",
+            output=cls.result_steering_distributed,
+        )
+        cls.runModule(
+            "r.mapcalc",
+            expression="ndvi_2002 = double(lsat7_2002_40@PERMANENT - lsat7_2002_30@PERMANENT) / double(lsat7_2002_40@PERMANENT + lsat7_2002_30@PERMANENT)",
+        )
+        cls.runModule(
+            "r.mapcalc",
+            expression="ndvi_1987 = double(lsat5_1987_40@landsat - lsat5_1987_30@landsat) / double(lsat5_1987_40@landsat + lsat5_1987_30@landsat)",
+        )
+        cls.runModule(
+            "r.mapcalc",
+            expression="urban_1987 = if(ndvi_1987 <= 0.1 && isnull(lakes), 1, if(isnull(lakes), 0, null()))",
+        )
+        cls.runModule(
+            "r.mapcalc",
+            expression="urban_2002 = if(ndvi_2002 <= 0.1 && isnull(lakes), 1, if(isnull(lakes), 0, null()))",
+        )
+        cls.runModule("r.slope.aspect", elevation="elevation", slope="slope")
+        cls.runModule("r.grow.distance", input="lakes", distance="lakes_dist")
+        cls.runModule("r.mapcalc", expression="lakes_dist_km = lakes_dist/1000.")
+        cls.runModule("v.to.rast", input="streets_wake", output="streets", use="val")
+        cls.runModule("r.grow.distance", input="streets", distance="streets_dist")
+        cls.runModule("r.mapcalc", expression="streets_dist_km = streets_dist/1000.")
+        cls.runModule(
+            "r.futures.devpressure",
+            input="urban_2002",
+            output="devpressure",
+            method="gravity",
+            size=15,
+            flags="n",
+        )
+        cls.runModule(
+            "r.watershed",
+            elevation="elevation",
+            drainage="drainage",
+            stream="streams",
+            threshold=1000,
+        )
+        cls.runModule(
+            "r.watershed", elevation="elevation", basin="basin", threshold=5000
+        )
+        cls.runModule("r.null", map="basin", null=1000)
+        cls.runModule(
+            "r.stream.distance",
+            stream_rast="streams",
+            direction="drainage",
+            elevation="elevation",
+            method="downstream",
+            difference="HAND",
+        )
+        cls.runModule("r.grow.distance", input="HAND", value="HAND_filled")
+        cls.runModule(
+            "r.lake",
+            elevation="HAND_filled",
+            water_level=8,
+            lake="flood_1",
+            seed="streams",
+        )
+        cls.runModule(
+            "r.lake",
+            elevation="HAND_filled",
+            water_level=4,
+            lake="flood_2",
+            seed="streams",
+        )
+        cls.runModule(
+            "r.lake",
+            elevation="HAND_filled",
+            water_level=2,
+            lake="flood_3",
+            seed="streams",
+        )
+        cls.runModule(
+            "r.mapcalc",
+            expression="flood_1 = if (flood_1, 0.05, null())",
+            overwrite=True,
+        )
+        cls.runModule(
+            "r.mapcalc",
+            expression="flood_2 = if (flood_2, 0.1, null())",
+            overwrite=True,
+        )
+        cls.runModule(
+            "r.mapcalc",
+            expression="flood_3 = if (flood_3, 0.2, null())",
+            overwrite=True,
+        )
+        cls.runModule(
+            "r.patch", input="flood_3,flood_2,flood_1", output="flood_probability"
+        )
+        cls.runModule("r.null", map="flood_probability", null=0)
+        cls.runModule(
+            "r.mapcalc", expression="acapacity = rand(-100, 100) * 0.01", seed=1
+        )
+        cls.runModule(
+            "r.mapcalc", expression="urban_2002_steered = if (urban_2002 == 0, -1, 0)"
+        )
+        cls.runModule(
+            "r.mapcalc", expression="urban_2002_flooded_steered = if (urban_2002 == 0, -1000, 0)"
+        )
+        rules_1 = """27511 = 27511
+                    27513 = 27513
+                    27518 = 27518
+                    27539 = 27539
+                    27606 = 27606
+                    27607 = 27607"""
+        gs.write_command(
+            "r.reclass",
+            input="zipcodes",
+            output="subregions_1",
+            rules="-",
+            stdin=rules_1,
+        )
+        rules_2 = """27529 = 27529
+                    27601 =27601
+                    27603 = 27603
+                    27604 = 27604
+                    27605 = 27605
+                    27608 = 27608
+                    27610 = 27610"""
+        gs.write_command(
+            "r.reclass",
+            input="zipcodes",
+            output="subregions_2",
+            rules="-",
+            stdin=rules_2,
+        )
+
+    @classmethod
+    def tearDownClass(cls):
+        cls.runModule(
+            "g.remove",
+            flags="f",
+            type="raster",
+            name=[
+                "slope",
+                "lakes_dist",
+                "lakes_dist_km",
+                "streets",
+                "streets_dist",
+                "streets_dist_km",
+                "devpressure",
+                "ndvi_2002",
+                "ndvi_1987",
+                "urban_1987",
+                "urban_2002",
+                "drainage",
+                "HAND",
+                "HAND_filled",
+                "streams",
+                "flood_1",
+                "flood_2",
+                "flood_3",
+                "flood_probability",
+                "acapacity",
+                "basin",
+                "urban_2002_steered",
+                "urban_2002_flooded_steered",
+                "subregions_1",
+                "subregions_2",
+                cls.result,
+                cls.result_steering,
+                cls.result_steering_distributed,
+            ],
+        )
+        cls.del_temp_region()
+        gs.try_remove(cls.output_demand)
+        gs.try_remove(cls.output_demand_1)
+        gs.try_remove(cls.output_demand_2)
+
+    def tearDown(self):
+        self.runModule(
+            "g.remove",
+            flags="f",
+            type="raster",
+            name=[
+                self.output,
+                self.output_1,
+                self.output_2,
+                self.output_development_pressure,
+                self.output_development_pressure_1,
+                self.output_development_pressure_2,
+            ],
+        )
+
+
+    def test_pga_run_simple_steering(self):
+        """Test if results is in expected limits"""
+        for i in range(0, 7):
+            self.assertModule(
+                "r.futures.simulation",
+                developed="urban_2002_steered" if i == 0 else self.output,
+                development_pressure="devpressure"
+                if i == 0
+                else self.output_development_pressure,
+                compactness_mean=0.4,
+                compactness_range=0.05,
+                discount_factor=0.1,
+                patch_sizes="data/patches.csv",
+                predictors=["slope", "lakes_dist_km", "streets_dist_km"],
+                n_dev_neighbourhood=15,
+                devpot_params="data/potential.csv",
+                random_seed=1,
+                num_neighbors=4,
+                seed_search="random",
+                development_pressure_approach="gravity",
+                gamma=1.5,
+                scaling_factor=1,
+                subregions="zipcodes",
+                demand="data/demand.csv" if i == 0 else self.output_demand,
+                output=self.output,
+                flags="r",
+                output_development_pressure=self.output_development_pressure,
+                output_demand=self.output_demand,
+                num_steps=1,
+                steering_step=i,
+                overwrite=True,
+            )
+        differences = compare_base_stats(self.result, self.output)
+        for difference in differences:
+            self.assertLess(difference, 5, "Cell count differs too much.")
+        self.assertRastersNoDifference(
+            actual=self.output, reference=self.result_steering, precision=1e-6
+        )
+
+    def test_pga_run_distributed_steering(self):
+        """Test if results is in expected limits"""
+        for i in range(0, 7):
+            self.assertModule(
+                "r.futures.simulation",
+                developed="urban_2002_steered" if i == 0 else self.output_1,
+                development_pressure="devpressure"
+                if i == 0
+                else self.output_development_pressure_1,
+                compactness_mean=0.4,
+                compactness_range=0.05,
+                discount_factor=0.1,
+                patch_sizes="data/patches.csv",
+                predictors=["slope", "lakes_dist_km", "streets_dist_km"],
+                n_dev_neighbourhood=15,
+                devpot_params="data/potential.csv",
+                random_seed=1,
+                num_neighbors=4,
+                seed_search="random",
+                development_pressure_approach="gravity",
+                gamma=1.5,
+                scaling_factor=1,
+                subregions="subregions_1",
+                demand="data/demand.csv" if i == 0 else self.output_demand_1,
+                output=self.output_1,
+                flags="r",
+                output_development_pressure=self.output_development_pressure_1,
+                output_demand=self.output_demand_1,
+                num_steps=1,
+                steering_step=i,
+                overwrite=True,
+            )
+            self.assertModule(
+                "r.futures.simulation",
+                developed="urban_2002_steered" if i == 0 else self.output_2,
+                development_pressure="devpressure"
+                if i == 0
+                else self.output_development_pressure_2,
+                compactness_mean=0.4,
+                compactness_range=0.05,
+                discount_factor=0.1,
+                patch_sizes="data/patches.csv",
+                predictors=["slope", "lakes_dist_km", "streets_dist_km"],
+                n_dev_neighbourhood=15,
+                devpot_params="data/potential.csv",
+                random_seed=1,
+                num_neighbors=4,
+                seed_search="random",
+                development_pressure_approach="gravity",
+                gamma=1.5,
+                scaling_factor=1,
+                subregions="subregions_2",
+                demand="data/demand.csv" if i == 0 else self.output_demand_2,
+                output=self.output_2,
+                flags="r",
+                output_development_pressure=self.output_development_pressure_2,
+                output_demand=self.output_demand_2,
+                num_steps=1,
+                steering_step=i,
+                overwrite=True,
+            )
+        self.runModule(
+            "r.patch", input=[self.output_1, self.output_2], output=self.output
+        )
+        differences = compare_base_stats(self.result, self.output)
+        for difference in differences:
+            self.assertLess(difference, 6, "Cell count differs too much.")
+        self.assertRastersNoDifference(
+            actual=self.output,
+            reference=self.result_steering_distributed,
+            precision=1e-6,
+        )
+
+    def test_pga_run_flooding_steering(self):
+        """Test if results is in expected limits"""
+        for i in range(0, 7):
+            self.assertModule(
+                "r.futures.simulation",
+                developed="urban_2002_flooded_steered" if i == 0 else self.output,
+                development_pressure="devpressure"
+                if i == 0
+                else self.output_development_pressure,
+                compactness_mean=0.4,
+                compactness_range=0.05,
+                discount_factor=0.1,
+                patch_sizes="data/patches.csv",
+                predictors=["slope", "lakes_dist_km", "streets_dist_km"],
+                n_dev_neighbourhood=15,
+                devpot_params="data/potential.csv",
+                random_seed=1,
+                num_neighbors=4,
+                seed_search="random",
+                development_pressure_approach="gravity",
+                gamma=1.5,
+                scaling_factor=1,
+                subregions="zipcodes",
+                demand="data/demand.csv" if i == 0 else self.output_demand,
+                output=self.output,
+                flags="r",
+                output_development_pressure=self.output_development_pressure,
+                output_demand=self.output_demand,
+                redistribution_matrix="data/matrix.csv",
+                flood_maps_file="data/flood_depth_input_steered.csv",
+                adaptive_capacity="acapacity",
+                huc="basin",
+                depth_damage_functions="data/damage_curves.csv",
+                ddf_subregions="zipcodes",
+                population_demand="data/population_demand.csv",
+                response_func=[0.25, 0.5, 0.25, 0.5],
+                response_stddev=0.1,
+                num_steps=1,
+                steering_step=i,
+                overwrite=True,
+            )
+
+    def test_pga_run_flooding_distributed_steering(self):
+        """Test if results is in expected limits"""
+        for i in range(0, 7):
+            self.assertModule(
+                "r.futures.simulation",
+                developed="urban_2002_flooded_steered" if i == 0 else self.output_1,
+                development_pressure="devpressure"
+                if i == 0
+                else self.output_development_pressure_1,
+                compactness_mean=0.4,
+                compactness_range=0.05,
+                discount_factor=0.1,
+                patch_sizes="data/patches.csv",
+                predictors=["slope", "lakes_dist_km", "streets_dist_km"],
+                n_dev_neighbourhood=15,
+                devpot_params="data/potential.csv",
+                random_seed=1,
+                num_neighbors=4,
+                seed_search="random",
+                development_pressure_approach="gravity",
+                gamma=1.5,
+                scaling_factor=1,
+                subregions="subregions_1",
+                demand="data/demand.csv" if i == 0 else self.output_demand_1,
+                output=self.output_1,
+                flags="r",
+                output_development_pressure=self.output_development_pressure_1,
+                output_demand=self.output_demand_1,
+                redistribution_matrix="data/matrix.csv",
+                flood_maps_file="data/flood_depth_input_steered.csv",
+                adaptive_capacity="acapacity",
+                huc="basin",
+                depth_damage_functions="data/damage_curves.csv",
+                ddf_subregions="zipcodes",
+                population_demand="data/population_demand.csv",
+                response_func=[0.25, 0.5, 0.25, 0.5],
+                response_stddev=0.1,
+                redistribution_external_output=self.redistribution_1,
+                num_steps=1,
+                steering_step=i,
+                overwrite=True,
+            )
+            self.assertModule(
+                "r.futures.simulation",
+                developed="urban_2002_flooded_steered" if i == 0 else self.output_2,
+                development_pressure="devpressure"
+                if i == 0
+                else self.output_development_pressure_2,
+                compactness_mean=0.4,
+                compactness_range=0.05,
+                discount_factor=0.1,
+                patch_sizes="data/patches.csv",
+                predictors=["slope", "lakes_dist_km", "streets_dist_km"],
+                n_dev_neighbourhood=15,
+                devpot_params="data/potential.csv",
+                random_seed=1,
+                num_neighbors=4,
+                seed_search="random",
+                development_pressure_approach="gravity",
+                gamma=1.5,
+                scaling_factor=1,
+                subregions="subregions_2",
+                demand="data/demand.csv" if i == 0 else self.output_demand_2,
+                output=self.output_2,
+                flags="r",
+                output_development_pressure=self.output_development_pressure_2,
+                output_demand=self.output_demand_2,
+                redistribution_matrix="data/matrix.csv",
+                flood_maps_file="data/flood_depth_input_steered.csv",
+                adaptive_capacity="acapacity",
+                huc="basin",
+                depth_damage_functions="data/damage_curves.csv",
+                ddf_subregions="zipcodes",
+                population_demand="data/population_demand.csv",
+                response_func=[0.25, 0.5, 0.25, 0.5],
+                response_stddev=0.1,
+                redistribution_external_output=self.redistribution_2,
+                num_steps=1,
+                steering_step=i,
+                overwrite=True,
+            )
+            if i < 6:
+                part_1_df = pd.read_csv(f"{self.redistribution_1}_{i + 1}.csv", index_col=0)
+                part_2_df = pd.read_csv(f"{self.redistribution_2}_{i + 1}.csv", index_col=0)
+                migrants_row = (part_1_df + part_2_df).sum()
+                demand_1_df = pd.read_csv(self.output_demand_1, index_col=0)
+                demand_2_df = pd.read_csv(self.output_demand_2, index_col=0)
+                # add all migrants to demand for next year
+                demand_1_df.iloc[i + 1] += migrants_row
+                demand_2_df.iloc[i + 1] += migrants_row
+                demand_1_df.round(0).astype(int).to_csv(self.output_demand_1)
+                demand_2_df.round(0).astype(int).to_csv(self.output_demand_2)
+
+        self.runModule(
+            "r.patch", input=[self.output_1, self.output_2], output=self.output
+        )
+
+
+if __name__ == "__main__":
+    test()