9, Final exercise TNO

Extending the RuimteScanner at TNO

Background

In modules 0 to 3 you learned the GeoDMS language: items and namespaces, units and domains, calculations across domains, classification and visualisation, reading and writing data, and meta scripting with templates and for_each. In this closing exercise, you apply all of that to a working model rather than a clean sheet.

The model is a simplified version of the RuimteScanner (Land Use Scanner, RS 1.0). It allocates a fixed set of land-use classes to grid cells. Three ingredients drive the allocation:

Ingredient What it is Where it lives in the config
Claims The regional demand (in hectares) per land-use class, per allocation region /Scenario_Components/Regional_Demand/<Scenario>/per_LU_type/<lu-class>
Suitability A score per grid cell per land-use class, expressing how fit a cell is for that use /Scenario_Components/Local_Suitability/<Scenario>
Allocation region The spatial unit at which claims are defined and balanced /ModelParameters/RegioUnit_str

The allocation engine then assigns each cell to the class that, given the claims and the suitability, produces the most plausible spatial pattern.

This exercise mirrors three concrete extensions that TNO wants to make to its own RS 1.0 configuration:

TNO goal Task in this exercise
Add more suitability layers Task 2
Define claims at a finer regional level (province, COROP, municipality) instead of Randstad, Intermediate and Periphery Task 3
Split land-use further, for example, economic sectors into sub-classes Task 4

By the end, you will have added a suitability factor, a finer allocation region, and a new land-use class to a running model, and you will have seen the effect of each change on the allocation result.

Before you start

You need a working GeoDMS GUI (module 0) and the exercise configuration.

  1. Download the TNO branch of the LUS demo here.
  2. Unzip it to your local project directory and open the main configuration in the GeoDMS GUI.
  3. The configuration runs out of the box with the original three-zone settings.

A note on how to work: extend in small steps and check the result of each step before moving on. After every change, recalculate and open a result item. If something turns red, the detail pages tell you what is missing. Build up, do not rewrite.

Getting oriented (no scripting yet)

Open the configuration and explore before you change anything.

  • Open the land-use classification Classifications/ggModel. How many classes are there? This is the list of land-use types the model can allocate.
  • Open the Suitability container. For one class, follow how its suitability score is built up from input layers. Notice that each contributing factor is a value on the grid that is normalised and then combined.
  • Open the Claims container. At which regional level are the claims currently defined? Open the claim table in a table view.
  • Find the model parameter RegioUnit. What does it point to now?
  • Open the final allocation result (/Final_Results/A1_GE_Discr) in a map view. This is the picture you will be changing.

Task 1: Reproduce the baseline (verify your setup)

Before extending anything, confirm that you can generate results with the current settings.

  1. Run the allocation with the settings as delivered (/Final_Results/A1_GE_Discr).
  2. Open the allocation result in a map view.
  3. Note two things you can compare against later: the overall spatial pattern and one number, for example, the allocated area of a single class.

This is your reference. Every later task changes the result, and you want to be able to see what changed.

Checkpoint: You have a result map and at least one reference number written down.

Task 2: Add a new suitability factor (Modules 1c, 1d, 2)

This is the most contained extension and a good warm-up. You add one new input layer and let it influence where a land-use class is allocated.

The supplied configuration contains an extra input layer in data that is not yet used, for example, a sand thickness layer. You will normalise it and add it to the suitability of one class.

  1. Add data. Open the new layer in the data/current_situation/thematic_maps/naturelandscape and view it on the grid. Make sure it is defined on the same domain as the other suitability inputs. If it is a polygon or point source, rasterise it to the grid first (module 2).
  2. Normalise it. A suitability contribution is a comparable score, not a raw value. Convert the raw layer to a score on a fixed range, for example, 0 to 1. Think about direction. If a high raw value means less suitable (for example, distance to a desired feature), invert the score with 1f - factor_score. And check your units: open your factor score in the Detail Pages and look at its values unit and metric. The raw layer carries the unit of whatever it measures; the score you feed into the suitability should not. If the metric still shows, say, metres, step 3 would be adding metres to scores, and this factor would drown out, or vanish next to, the others.
  3. Wire it into suitability. Open the suitability definition of the class you want to influence (for example, the residential or the working class) and add your normalised factor as an extra weighted term.
  4. Run and compare. Recalculate, open the allocation result, and compare it to your Task 1 reference. Where did the allocation of this class move towards or away from?

Reflection: change the weight up and down and rerun. How sensitive is the pattern to this one factor? This is the same tuning judgement you will need for any real suitability layer.

Checkpoint: the new factor appears in the suitability map of your chosen class, and the allocation result has visibly shifted.

Task 3: Add a finer allocation region (Modules 1c, 2)

TNO wants claims at the province, COROP or municipality level instead of the three zones. Here, you switch the allocation region to a finer CBS area classification and supply claims at that level. Download a geopackage of the regions here.

  1. Add data to /Geography/RegionalUnits. Confirm the finer region (for example, province) is read in and rasterised to the grid, so that every cell has a region. Open it on the map.

  2. Adjust the model parameter. Point RegioUnit at the finer region. Recalculate and follow where the model now expects claims. It will report missing claims at the new level, which is the next step.

  3. Provide claims at the new level. You have two options:

    Option A, author them. Make a claim table indexed by the new region and fill it with plausible demands, just like the current Excel file. Option B, use the provided helper. The demo includes a helper that derives finer-region claims from the original three-zone claims by splitting each zone claim over its regions in proportion to the current area of that class per region. The totals per class are preserved, so the model stays balanced exactly as before, only at a finer spatial resolution. The core of it looks like this, per class:

// current area of a class per region
attribute<claimHa> Current_Provincie (Provincie) := ='sum(Current_situation/current_landuse/'+lu_name+', provincie/per_domain)';

// share of the region within its zone, via the zone relation
attribute<float64> Share_Provincie (Provincie) := Current_Provincie / Current_RandstadIntermediatePeripheral[Provincie/Landsdelen_rel];

// finer claim = original zone claim of that region times its share
attribute<claimHa> Claim_Provincie (Provincie) := Claim_RandstadIntermediatePeripheral[Provincie/Landsdelen_rel] * Share_Provincie;

// claims are additional, so add the current area to it.
attribute<claimHa> Total_Provincie (Provincie) := Current_Provincie + Claim_Provincie;

For the demo, option B gives you a working result without hand-authoring a single number.

Answer for B 
Template Claims_per_LU_T
{
	parameter<string>    lu_name;
	attribute<RegioUnit> Regios (domain);
	///
	

	attribute<claimHa>      Current_RandstadIntermediatePeripheral (RandstadIntermediatePeripheral) := ='sum(Current_situation/current_landuse/'+lu_name+', RandstadIntermediatePeripheral/per_domain)';
	attribute<claimHa>      Claim_RandstadIntermediatePeripheral   (RandstadIntermediatePeripheral) := ='rjoin(id(RandstadIntermediatePeripheral)+1,Read_Excel/'+lu_name+'/RstIntPer[uint32], Read_Excel/'+lu_name+'/'+scenario_name_short+'[claimHa])'; //claims in excel hebben regio=0 voor exterior. RegioUnit heeft geen exterior, dus +1
	attribute<claimHa>      Total_RandstadIntermediatePeripheral   (RandstadIntermediatePeripheral) := Current_RandstadIntermediatePeripheral + Claim_RandstadIntermediatePeripheral;
	
	// hack om landsdeel claims te converteren naar Provinces
	attribute<claimHa>      Current_Province                                             (Province) := ='sum(Current_situation/current_landuse/'+lu_name+', Province/per_domain)';
	attribute<float64>      Share_Province                                               (Province) := Current_Province / Current_RandstadIntermediatePeripheral[Province/Landsdelen_rel];
	attribute<claimHa>      Claim_Province                                               (Province) := Claim_RandstadIntermediatePeripheral[Province/Landsdelen_rel] * Share_Province;
	attribute<claimHa>      Total_Province                                               (Province) := Current_Province + Claim_Province;
	
	attribute<claimHa>      Claim                                                       (RegioUnit) := ='Claim_'+ModelParameters/RegioUnit_str;
	attribute<claimHa>      Total                                                       (RegioUnit) := ='Total_'+ModelParameters/RegioUnit_str;
	attribute<OperatorType> Operator                                                    (RegioUnit) := const(1, RegioUnit, OperatorType);
	attribute<EUR_M2>       logAj                                                       (RegioUnit) := const(0, RegioUnit, Eur_m2);
}


----

unit<UInt32> Province := CBS_Gebiedsindelingen/Provincie_gegeneraliseerd
, Descr = "CBS Provincie indeling"
{
	attribute<rdc_meter>      geometry                                 (poly) := CBS_Gebiedsindelingen/Provincie_gegeneraliseerd/geometry;
	attribute<String>         Label                                           := statnaam;
	attribute<String>         name                                            := replace(AsItemName(statnaam), 'Fryslân', 'Friesland');
	attribute<.>              Per_rdc_100m                         (rdc_100m) := poly2grid(geometry, rdc_100m);
	attribute<.>              per_domain                             (domain) := poly2grid(geometry, domain);
	
	attribute<RandstadIntermediatePeripheral>     Landsdelen_rel := switch(
													case(id(.) == v/Groningen || id(.) == v/Drenthe || id(.) == v/Zeeland || id(.) == v/Limburg || id(.) == v/Friesland || id(.) == v/Overijssel, RandstadIntermediatePeripheral/v/Peripheral)
													, case(id(.) == v/Gelderland  || id(.) == v/Noord_Brabant || id(.) == v/Flevoland, RandstadIntermediatePeripheral/v/Intermediate)
													, case(id(.) == v/Noord_Holland || id(.) == v/Zuid_Holland || id(.) == v/Utrecht, RandstadIntermediatePeripheral/v/Randstad)
													, null_u
												);
	container V := for_each_nedv(name, string(id(.))+'[..]', void, .);
}

unit<UInt32> COROP     := CBS_Gebiedsindelingen/coropgebied_gegeneraliseerd
, Descr = "CBS COROP indeling"
{
	attribute<rdc_meter>      geometry                                 (poly) := CBS_Gebiedsindelingen/coropgebied_gegeneraliseerd/geometry[rdc_meter];
	attribute<String>         Label                                           := statnaam;
	attribute<String>         name                                            := AsItemName(statnaam);
	attribute<.>              Per_rdc_100m                         (rdc_100m) := poly2grid(geometry, rdc_100m);
	attribute<.>              per_domain                             (domain) := poly2grid(geometry, domain);
	
	attribute<RandstadIntermediatePeripheral>     Landsdelen_rel := point_in_polygon(centroid_or_mid(geometry), RandstadIntermediatePeripheral/geometry);
}
  1. Run and compare. Recalculate and open the result. The regional pattern of the claims is now finer. Compare to Task 1: the same total demand is now distributed across more regions.

Reflection: try a different CBS level (COROP or municipality) by changing the RegioUnit and rerunning. What happens to the result as the regions get smaller? Where do you start to see the limits of the available data?

Checkpoint: the allocation runs at the finer region level, and the result reflects the finer claims.

Task 4: Add a new land-use sector (Modules 1c, 1d, 3)

This is the most involved extension. TNO wants to split land-use further, for example, an economic sector into subclasses. Here, you add one new class to the model end-to-end. The same recipe works for any split.

A new class touches four parts of the model: the classification, its start state, its suitability, and its claim. Do them in this order and check after each.

  1. Extend the classification. Add the new class to Classifications/ggModel: a class code, a name, and a BrushColor. Recalculate and confirm the class appears in the classification table. In this example, add an extra class called ‘services’, which splits the BBG-service classes from residential into its own class.

  2. Supply its start state. The model needs to know where the class currently exists on the grid. In this demo, the start-state maps are provided here, so add the supplied maps to the folder: %projdir%/data/Current_situation/bbg_2022. Open it on the map and check that it lands in plausible places.

  3. Add its suitability. Define a suitability score for the new class, following the pattern of an existing class. Reuse the factors that make sense for it and adjust the weights.

  4. Add its claim. Add a claim for the new class at the active region level. If you are using the helper from Task 3, add the new class to that split so it gets a finer-region claim too.

  5. Run and verify. Recalculate and open the result. Confirm the new class appears in the allocation, in roughly the right places, and that its allocated total is in line with the claim you gave it.

Reflection and a caveat. RS 1.0 allocates one class per cell, so a split is only as fine as your start-state and suitability data allow. To split, say, an industry into sub-sectors by economic activity (SBI), you need data that distinguishes those sub-sectors in space, for example, an establishment or terrain typology. Where such data does not exist, the split becomes a proxy rather than a true sub-allocation. This is a real boundary of the 1.0 approach and a good thing to discuss as a group.

Checkpoint: the new class is allocated and its total matches its claim.

Task 5 (open): Apply it to your own question

Pick one extension you would actually want in your own TNO configuration.

  • Which region level do you need, and is the data available at that level?
  • Which class would you split, and what data distinguishes the sub-classes in space?
  • Which suitability factor is missing for your case, and where would you get it?

Sketch the four touch points (classification, start state, suitability, claim, or for a region: data, parameter, claims) for your own case. This is the bridge from the demo to your own work.

What you have practised

Task Topic RuimteScanner part
1 Running the model, reading results Allocation result
2 Normalising and combining input layers Suitability
3 Switching the allocation region, deriving claims Allocation region and claims
4 Adding a land-use class end to end Classification, start state, suitability, claims
5 Translating the recipe to your own case All of the above

Go to previous module: 8, Final exercise Deltares

Go to next module: Module 10, Starting your own project