The NetLogo GIS extension: a tutorial

I have previously written about the benefits of extensions to NetLogo. One of the powerful built-in extensions to NetLogo is the Geographic Information Systems (GIS) extension, which allows the import of vector and raster GIS data into NetLogo simulation models. This is particularly helpful in building simulations of real populations (as well as allowing the import of such things as terrain data).

In this tutorial, we discuss a modified version of Uri Wilensky’s “GIS General Examples” model, which is bundled with NetLogo in the “Code Examples” section of the Models Library. This tutorial example requires a copy of the data folder associated with that model (which can be found in the “Models” folder of the NetLogo installation directory). In particular, this tutorial example imports country shape data and (slightly old) country population data, so as to be able to compute population densities, like so:

The data behind this includes the following table:

CNTRY_NAME	POP_CNTRY	SQKM	…
Afghanistan	17250390	641358.44	…
Albania	3416945	28798	…
Algeria	27459230	2323510.25	…
…	…	…	…

As well as colouring countries by population density, we will be able to construct representative models of the world population, where each NetLogo agent represents 10 million people in the real world, and these agents are located in the appropriate country on a world map. For example, Australia’s population of roughly 20 million people is represented by 2 agents, while Canada’s population of roughly 30 million people is represented by 3 agents. India and China are a little more crowded:

Constructing agent populations in this way allows agents to have country-dependent behaviour. The same approach, with different GIS datasets, allows models of individual countries, states, or cities, with populations determined by GIS data on, for example, counties or suburbs.

The NetLogo code

Loading the GIS data is much like Uri Wilensky’s “GIS General Examples” model (specific details will be given later), while the following procedure sets up the population density data in the patches:

to setup-gis
  show "Loading patches..."
  gis:apply-coverage countries-dataset "POP_CNTRY" population
  gis:apply-coverage countries-dataset "SQKM" area
  gis:apply-coverage countries-dataset "CNTRY_NAME" country-name

  ask patches [
    ifelse (area > 0 and population > 0)
      [ ; Colour patch according to population density
        set population-density (population / area)
        set pcolor (scale-color red population-density 400 0) ]
    
      [ ; Colour patch with no population blue
        set population-density 0
        set pcolor blue ]
  ]
end

Here the calls to gis:apply-coverage copy values from the GIS country dataset to the named patch variables. The ask which follows calculates the population density of each patch, and colours it appropriately.

How should we create the agents? The simplest answer is to use those calculated population densities, combined with the fact that the average patch has an area of 510,000,000 / N square kilometres, where N is the number of patches (ignoring the substantial distortion due to the map projection). To take India as an example, the population density is 283.7 people per square kilometre, which works out to 13,671,000 people per patch, which will be represented by 1.3671 agents (this is in fact a slight overestimate, since the sum total of India’s patches have a slightly greater area than the country itself). We can effectively create 1.3671 agents per patch by always creating one agent, and creating an additional one with probability 0.3671. The sprout-<breeds> primitive is used to create the agents:

to create-random
  ; Delete old population
  ask persons [ die ]

  ask patches [
    if (population-density > 0) [
      ; How many people to create? (N.B. this is an overestimate, since country does not occupy entire patch)
      let num-people (population-density * patch-area-sqkm / 10000000)
      
      ; Create whole numbers directly
      sprout-persons (floor num-people) [ agent-setup country-name ]
      
      ; Create fractions probabilistically
      let fractional-part (num-people - (floor num-people))
      if (fractional-part > random-float 1) [ sprout-persons 1 [ agent-setup country-name ] ]
    ]
  ]
  
  show (word "Randomly created " (count persons) " people")
  tick
end

to agent-setup [ ctry ] ; set up agent parameters, including the country of the agent
  set shape "person"
  set size 3
  set color black
  set agent-country-name ctry
end

Results

This code gives a different result each time, with between 540 and 640 agents created. The black line on the histogram below shows the average value of 584, while the dotted line shows the actual world population (some years ago) of 5.695 billion (corresponding to 570 agents). Notice that over-estimating the number of agents per grid square has given a population which is too large. This is in spite of the fact that populations of countries occupying zero grid squares don’t get counted at all. Examples of such small countries include Singapore (population 2.8 million), Rwanda (7.9 million), and Burundi (6.0 million). The other difficulty is that, because of the substantial distortion due to the map projection, there are not enough agents in countries near the equator (like Brazil), and too many in countries near the poles (like Russia). The dashed blue line shows the result of a better (and deterministic) method, which we discuss next.

A better way

For the deterministic alternative, we also use the NetLogo tables extension, which provides tables of key-value mappings (with put, get, and has-key?). We use two tables: patches-to-go, which maps a country name to a number of patches, and people-to-go, which maps a country name to a number of agents. The second of these tables is created using a loop over GIS attributes (3):

to create-deterministic
  ; Reset tables (1)
  set patches-to-go table:make
  set people-to-go table:make
  set desired-people 0
  
  ; Delete old population (2)
  ask persons [ die ]
  
  ; Calculate number of people desired (3)
  ; N.B. Countries with no patches will be ignored, even if population > 0
  foreach gis:feature-list-of countries-dataset [
    let ctry gis:property-value ? "CNTRY_NAME"
    let pop gis:property-value ? "POP_CNTRY"
    let desired (round (pop / 10000000))
    table:put people-to-go ctry desired
  ]

This sets up the people-to-go table to hold the number of agents that must be created for each country (the question marks refer to the “vector feature” objects that the loop is iterating over). The number of patches for each country, on the other hand, is trickier. It is tempting to use multiple calls to count patches with ..., but that would be slow. The following alternative iterates once over all patches, adding one to the appropriate table entry in patches-to-go each time (and also accumulating the total desired-people to create, excluding countries with zero patches):

  ; Count patches for each country, efficiently, using a single loop (4)
  ask patches [
    if (population-density > 0) [
      
      ; Have we seen this country before?
      ifelse (table:has-key? patches-to-go country-name)
        [ ; YES -- record one extra patch for this country
          let n (table:get patches-to-go country-name)
          table:put patches-to-go country-name (n + 1) ]
        
        [ ; NO -- record first patch, and add population of country to desired total
          table:put patches-to-go country-name 1
          let desired (table:get people-to-go country-name)
          set desired-people (desired-people + desired) ]
    ]
  ]

The actual creation process looks at the ratio of the people-to-go and patches-to-go table entries, creating agents when this ratio becomes large enough. As the dashed blue line on the histogram above indicates, the number of agents created (558) is less than the 570 expected from the total world population. This is because of the exclusion of countries with zero patches, as well as the impact of rounding. Making the patches smaller alleviates this problem (with patches of a tenth the area, 563 agents are created), but slows down execution.

  ; Create desired number of people (5)
  ask patches [
    if (population-density > 0) [
      let n (table:get patches-to-go country-name)
      let people-needed (table:get people-to-go country-name)
      
      if (people-needed > 0) [
        ; More people needed for this country
        table:put patches-to-go country-name (n - 1)
        let num-people (round (people-needed / n))
        
        if (num-people > 0) [
          ; People needed on this patch
          table:put people-to-go country-name (people-needed - num-people)
          sprout-persons num-people [ agent-setup country-name ]
        ]
      ]
    ]
  ]
  
  ; Report number of agents created (6)
  show (word "Created " (count persons) " people / desired " desired-people)
  tick
end

Some basics

Finally, the following setup code loads the GIS data, much as in Uri Wilensky’s “GIS General Examples” model. As with many NetLogo features, there is considerable sophistication “under the hood,” making the GIS data import relatively easy. However, the gis:load-dataset procedure will not, unfortunately, load all SHP files.

One point to note is the use of the gis:set-world-envelope-ds primitive to fit the GIS data to the NetLogo world grid. The “DS” version of this primitive is necessary because having max-pxcor twice the size of max-pycor does not give the world the exact 2:1 aspect ratio which the dataset uses.

extensions [ gis table ]
globals [ projection countries-dataset patch-area-sqkm desired-people patches-to-go people-to-go ]
breed [ persons person ]
patches-own [ population country-name area population-density ]
persons-own [ agent-country-name ]

to setup
  ca
  set projection "WGS_84_Geographic"
  
  ; Note that setting the coordinate system here is optional
  gis:load-coordinate-system (word "data/" projection ".prj")
  
  ; Load the dataset
  set countries-dataset gis:load-dataset "data/countries.shp"
  
  ; Set the world envelope. DS = different x and y scales
  ; (fits to actual world dimensions -- world does not actually have a perfect 2:1 aspect ratio)
  gis:set-world-envelope-ds [-180 180 -90 90]

  ; Drawing country boundaries from a shapefile -- as per Uri Wilensky
  gis:set-drawing-color white
  gis:draw countries-dataset 1
  
  ; New code
  set patch-area-sqkm (510000000 / count patches)
  setup-gis
  reset-ticks
end

Going further

Since each agent contains a country name, we can extend this tutorial example with country-specific economic data or population-growth data, and so simulate evolution of the state of the planet using country-specific behaviour. For anyone wishing to experiment with doing so, the full NetLogo code for this tutorial can be downloaded from Modeling Commons. As an example, we can use NetLogo file-read operations to read in this list of countries and median ages into another table (data adapted from Wikipedia). The following additional code in agent-setup will then assign a moderately realistic age to each agent (assuming ages very roughly follow a uniform distribution, and ignoring the differences between mean and median):

  let med-age 28.4  ; use overall world median as default
  if (table:has-key? country-median-age ctry) [
    set med-age (table:get country-median-age ctry)
  ]
  set age random-float (2 * med-age)

There are also complexities associated with building models satisfying multiple constraints (e.g. combining the kind of spatial constraints discussed here with various demographic profiles). For information on that, I refer readers to the agent-based modelling literature (e.g. this book).

Readers can also experiment with loading other GIS datasets (although the GIS extension will not load all SHP files). For example, here is the result of loading Statistical Local Areas from the Australian Bureau of Statistics:

Advertisement