Revisiting Artificial Anasazi: a NetLogo tutorial (part 1)

I have previously uploaded some NetLogo tutorials, and for the next one I thought I would revisit the famous “Artificial Anasazi” model of Axtell et al. (2002). This model was ported to NetLogo by Marco Janssen in 2008 (see his 2009 JASSS paper and his model). I have rewritten the model from scratch for explanatory purposes (though very much inspired by Janssen’s work). It makes a good example of how NetLogo can assist land-use studies, both archaeological and present-day.

The model studies the population of Anasazi farmers living in Long house Valley (near Black Mesa in Arizona) between 800 and 1350 AD. The map below shows part of the valley, through which Highway 160 now runs (there are more maps here):

Much has been preserved by the dry Arizona climate. Fragments of pottery (like the bowl at the top of the page) are found routinely, and these well-preserved ruined Anasazi dwellings below are not far away. Extensive archaeological investigation has told us both the climate that existed between 800 and 1350 AD, and the population of the valley during that time.

For purposes of simplicity, my version of the model does not include the availability of drinking water and the locations of dwellings. I am only modelling the location of farms. Consequently, although I hope to accurately model population numbers, I do not hope to accurately model the locations at which people lived, as Axtell et al. did. In addition, this first half of the tutorial will only look at code for the patches in the model.

The model requires three data files in a folder called “adata” (zipped version here). One of these is a map of the valley, formatted as an ascii GIS raster file, which can be read by the NetLogo GIS extension. Such use of GIS datasets is common in land-use models. The following command will read in the raster file, and set a numerical value attribute of each patch.

to gis-map-load [ fname ]
  ifelse (file-exists? fname)
    [ let dataset gis:load-dataset fname
      gis:set-world-envelope (gis:envelope-of dataset)
      gis:apply-raster dataset value
    ]
    [ file-error fname ]
end

to file-error [ fname ]
  user-message (word "Cannot find the file \"" fname "\"")
end

These numerical value attributes encode different regions of the valley. The following two reporters (functions) transform numerical codes to names (which will be stored in a zone attribute), and transform names to colours (which will be stored in the standard patch colour attribute):

to-report patch-zone-name [ i ]
  ifelse (i = 0) [ report "General" ] [      ;; General valley floor
  ifelse (i = 10) [ report "North" ] [       ;; North valley floor
  ifelse (i = 15) [ report "North Dunes" ] [ ;; North valley dunes
  ifelse (i = 20) [ report "Mid" ] [         ;; Mid valley floor
  ifelse (i = 25) [ report "Mid Dunes" ] [   ;; Mid valley dunes
  ifelse (i = 30) [ report "Natural" ] [     ;; Natural (non-arable)
  ifelse (i = 40) [ report "Uplands" ] [     ;; Uplands (arable)
  ifelse (i = 50) [ report "Kinbiko" ] [     ;; Kinbiko Canyon
  ifelse (i = 60) [ report "Empty" ] [ 
  report "?" ] ] ] ] ] ] ] ] ]               ;; there should be no "?" patches
end

to-report patch-zone-color [ z ]
  ifelse (z = "General") [ report brown ] [
  ifelse (z = "North") [ report red ] [
  ifelse (z = "Mid") [ report gray ] [
  ifelse (z = "Natural") [ report yellow ] [
  ifelse (z = "Uplands") [ report orange ] [
  ifelse (z = "Kinbiko") [ report pink ] [ 
  ifelse (z = "North Dunes") [ report green ] [ 
  ifelse (z = "Mid Dunes") [ report green ] [ 
  ifelse (z = "Empty") [ report white ] [ 
  report magenta ] ] ] ] ] ] ] ] ]  ;; there should be no magenta patches
end

I have used a formatting style here which is suitable for highly nested ifelse structures in NetLogo (the NetLogo editor helps ensure the correct number of right brackets). The following utility code adds a label to a patch, and loading the GIS file followed by some labelling produces the map below:

to make-region-label [ x y txt clr ]
  ask patch x y [
    set plabel-color clr
    set plabel txt
  ]
end

to setup  ;; this code will be updated later in the tutorial
  clear-all
  gis-map-load "adata/Map.asc"
  ask patches [
    set zone (patch-zone-name value)
    if (zone = "?") [ user-message (word "Error in patch data load for " self) ]
    set pcolor (patch-zone-color zone)
  ]
  
  make-region-label 45 98  "Non-arable uplands" black
  make-region-label 28 47  "General valley" brown
  make-region-label 23 73  "North valley" red
  make-region-label 50 68  "Mid valley" gray
  make-region-label 25 84  "Dunes" green
  make-region-label 20 102 "Kinbiko Canyon" pink
  make-region-label 39 11  "Arable uplands" orange
  ...
end

Here is the colour-coded NetLogo map of the valley:

The colour-coding is important, because different parts of the valley have had different fertility levels over time. Fortunately, we have a good idea what those levels were. One of the data files contains estimates of the Palmer Drought Severity Index (PDSI) for each patch, over the 800–1350 period of the simulation. This file is formatted as a sequence of lists, one for each different zone of the valley, with the first element of each list being the zone name, and the rest of the list being PDSI numbers (it’s handy that NetLogo can read in an entire list in one gulp). The following code transforms the data file into a table called pdsi-table (obviously this also requires the NetLogo “table” extension).

to-report load-zone-file [ fname ]
  let tbl table:make
  ifelse (file-exists? fname)
    [ file-open fname
      while [ not file-at-end? ] [
        let lst file-read
        table:put tbl (item 0 lst) (but-first lst)  ;; split first item (zone name) from rest of list
      ]
      file-close
    ]
    [ file-error fname ]
  report tbl
end

to setup
  ...
  set pdsi-table load-zone-file "adata/ZoneAdjPDSI.txt"
  set start-year 800
  set year start-year
  ...
end

The data file used here was extracted from that used in the original model. Having read it, the following simple reporter (function) can then look up the PDSI value for a particular zone and year:

to-report get-pdsi [ z y ]
  ifelse (table:has-key? pdsi-table z)
    [ report item (y - start-year) (table:get pdsi-table z) ]
    [ report 0 ]
end

For each zone, the PDSI determines an estimated crop yield for the farmers. Although the Anasazi grew a range of crops (corn, beans, and squash – see coin above), the model expresses crop yields as if only corn was grown. In the original model, the relationship between PDSI and crop yield is expressed as a step function, but I have smoothed this out using a cubic polynomial interpolation (the red curve below, rather than the dashed step function):

The following code calculates the crop yields for each patch (in kilograms per hectare), using these cubic interpolations:

to-report thresholded-cubic [ x a b c d lo hi ]
  let res (x * x * x * a + x * x * b + x * c + d)
  if (res < lo) [ set res lo ]
  if (res > hi) [ set res hi]
  report res
end

to calculate-patch-yields  ;; calculate the crop yield for each patch based on the PDSI data, using a cubic interpolation of published numbers
  ask patches [
    let pdsi get-pdsi zone year
    let theoretical-yield 0

    if (zone = "North" or zone = "Kinbiko" or (zone = "Mid" and pxcor <= 34)) [
      set theoretical-yield (thresholded-cubic pdsi 4.4167 6.9456 49.583 823.48 617 1153)
    ]
    if (zone = "General" or (zone = "Mid" and pxcor > 34)) [
      set theoretical-yield (thresholded-cubic pdsi 3.65 5.7925 41.65 686.28 514 961)
    ]
    if (zone = "Uplands") [
      set theoretical-yield (thresholded-cubic pdsi 2.9333 4.6599 33.267 548.77 411 769)
    ]
    if (zone = "North Dunes" or zone = "Mid Dunes") [
      set theoretical-yield (thresholded-cubic pdsi 4.5833 7.1871 51.917 858.03 642 1201)
    ]

    set base-yield (theoretical-yield * ...)  ;; missing code to be supplied later
    ...
  ]
end

The base-yield attribute of each patch is based on the theoretical yield calculated in this way, but it is modified in two ways, using the lower two of the following sliders:

First, we introduce some spatial variability in crop yields, using a normally distributed patch-quality attribute:

to-report apply-variability [ var ]
  let res random-normal 1 var  ;; random variation around a mean of 1
  if (res < 0) [ set res 0 ]   ;; because random-normal can return negative values
  report res
end

to setup
  ...
  ask patches [
    set patch-quality (apply-variability crop-yield-variability)
  ]
  ...
end

Because the normal distribution (the “bell curve”) has infinitely long tails, the calls to random-normal can sometimes return negative or very large crop yield values. Excessively large values have little effect on the simulation, because the model incorporates a rule that stored corn more than two years old goes bad and is thrown away. However, negative values must be explicitly tested for (negative crop yields are likely to give some rather strange behaviour!).

The second adjustment is that the entire crop yield is discounted using the harvest-adjustment slider:

to calculate-patch-yields
  ask patches [
    ...
    set base-yield (theoretical-yield * patch-quality * harvest-adjustment) 
    ...
  ]
end

In the models of Axtell et al. and Janssen, low values of harvest-adjustment are used to avoid unrealistically high population levels. The histogram below shows a typical run, for the year 1200 AD (a good year) on the General valley floor. Here, the theoretical yield of 961 is discounted to a mean of 560 kg (dashed line), well below the 800 kg (dotted line) required to support a typical five-person household for a year:

The effect of the harvest-adjustment is thus to restrict the population of the General valley floor to 53 good farms (8.3% of the patches), but also to impose a rather unrealistic triangular crop-yield distribution on the patches being farmed (i.e. the right-hand tail of the “bell curve”). It seems to me better to explicitly provide a usable-farm-fraction slider that restricts farming to a subset of patches. This permits a more realistic harvest-adjustment of 0.9, corresponding to “wastage” of only 10% of the crop.

It turns out that 0.21 is a good value for usable-farm-fraction. But why would 79% of the patches be unavailable for farming? One reason is the factors I have excluded in this tutorial version of the model – water supply and housing. Some potential farms are just too far away from suitable housing sites. The Anasazi penchant for cliff dwellings also suggests that potential housing sites had to be defensible, and this also rules out some potential farms as being too far away from suitable housing. Additionally, some patches would be ruled out by various terrain factors, and there would also be social pressures that limited the population level. Furthermore, the patches in the model are a little less than a hectare in size, and so one-hectare farms will occupy a little more than one patch. The usable-farm-fraction slider compensates for all of these effects.

We package the various patch initialisation activities together as follows. This includes using usable-farm-fraction to select a set of patches which will have non-zero patch-quality (as usual, we obtain a certain probability by testing for random-float 1 being less than that probability). The code also includes a historical-base-yield attribute, which will hold a list of the three most recent base yields for a patch:

to initialise-patch
  set zone (patch-zone-name value)  ;; as described above
  if (zone = "?") [ user-message (word "Error in patch data load for " self) ]
  set pcolor (patch-zone-color zone)
  set being-farmed? false
  set historical-base-yield [ ]  ;; this will be expanded to a list of the last 3 yields
  
  ifelse (random-float 1 < usable-farm-fraction)
    [ set patch-quality (apply-variability crop-yield-variability) ]
    [ set patch-quality 0 ]
end

We call this initialisation in setup, and we update the historical-base-yield list in the calculate-patch-yields procedure. We can also average the list to smooth out temporal variation in crop yields, compensating for the corn-storage ability of the Anasazi. The code count patches with [ mean historical-base-yield >= 800 ] will count feasible farms, taking corn storage into account.

to setup
  ...
  ask patches [
    initialise-patch
  ]
  ...
end

to calculate-patch-yields
  ask patches [
    ...  ;; calculate theoretical-yield as described above
    set base-yield (theoretical-yield * patch-quality * harvest-adjustment)
    set historical-base-yield (fput base-yield historical-base-yield)  ;; put new base-yield on the front of the list
    if (length historical-base-yield > 3) [ set historical-base-yield (but-last historical-base-yield) ]  ;; truncate to 3 entries, if necessary
  ]
end

The main additional code in setup is creation of a list of historical population levels:

to setup
  clear-all  ;; as previously described
  gis-map-load "adata/Map.asc"
  set pdsi-table load-zone-file "adata/ZoneAdjPDSI.txt"
  ask patches [ initialise-patch ]
  set start-year 800
  set year start-year

  set historical-sequence
    [ 14 14 14 14 14 14 14 14 14 14 13 13 13 12 12 12 12 12 12 12 11 11 11 11 11 10 10 10 10 10 9 9 9 9 9 9 9 9 8 8 7 7 7 7 7 7 7 7 7 7
      28 29 29 29 29 29 29 29 29 29 29 30 30 31 31 31 31 32 32 32 32 33 33 33 33 33 37 37 37 37 37 39 39 39 40 40 40 40 41 41 41 42 42
      42 42 42 42 42 42 42 56 58 58 58 58 58 58 58 58 58 58 60 60 60 60 60 60 61 61 61 60 61 61 61 61 59 61 61 61 61 61 62 62 62 62 62
      62 62 62 62 61 63 63 63 63 63 63 63 63 63 66 67 67 67 67 67 67 67 67 67 66 67 67 67 67 67 67 66 66 66 66 69 69 69 69 65 68 68 68
      68 66 67 67 67 66 66 66 66 66 66 66 68 68 68 68 68 68 68 68 68 87 87 87 87 87 87 87 86 86 87 85 86 86 87 87 87 87 88 88 88 87 87
      87 87 87 88 90 90 90 90 89 91 92 92 95 95 95 95 97 97 94 94 94 94 95 95 95 95 95 95 134 139 139 139 139 139 139 142 142 142 142
      143 143 146 146 146 146 151 151 153 151 151 151 151 151 156 164 164 164 164 163 163 163 163 165 165 165 165 164 164 163 164 164
      164 161 162 162 162 162 162 161 164 164 165 165 166 166 166 167 166 166 166 166 159 159 160 160 158 159 160 160 161 162 162 162
      148 150 151 151 151 151 151 149 149 147 148 148 148 148 150 150 150 151 151 152 152 152 153 153 153 116 117 118 118 119 119 120
      121 122 126 124 124 124 126 127 127 129 131 131 133 132 134 134 136 137 133 138 138 139 139 138 140 140 142 142 142 143 143 144
      145 145 147 146 147 147 149 149 150 151 151 145 146 146 147 148 149 149 151 153 154 155 157 159 160 161 163 164 166 167 167 169
      169 171 171 173 170 173 176 176 178 178 180 182 184 184 185 188 189 190 192 191 193 192 194 194 194 195 196 199 200 192 193 194
      196 196 197 200 202 202 204 201 209 208 211 212 212 213 214 215 216 210 208 206 201 196 188 181 176 172 167 159 156 148 146 141
      120 118 114 106 103 95 90 88 83 74 70 68 60 58 56 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
      0 0 0 0 0 0 0 0 0 0 0 0 ]
  
  make-region-label 45 98  "Non-arable uplands" black  ;; as previously described
  make-region-label 28 47  "General valley" brown
  make-region-label 23 73  "North valley" red
  make-region-label 50 68  "Mid valley" gray
  make-region-label 25 84  "Dunes" green
  make-region-label 20 102 "Kinbiko Canyon" pink
  make-region-label 39 11  "Arable uplands" orange
 
  reset-ticks
  tick-advance start-year  ;; start counting from the year 800
end

The code for a step just plots crop yields (smoothed, using the historical-base-yield list) and calculates new crop yields:

to go  ;; code for one step
  calculate-patch-yields
  plot-interesting-data
  if (year = 1350) [ stop ]
  set year year + 1
  tick
end

to plot-interesting-data  ;; plotting code
  set-current-plot "Population (households)"
  set-current-plot-pen "Potential"
  plotxy year (count patches with [ mean historical-base-yield >= 800 ])
  set-current-plot-pen "Historical"
  plotxy year (item (year - 800) historical-sequence)  ;; index list starting at 0
end

Plotting (see graph below) shows a good fit between predicted carrying capacity of the valley and actual population levels. Even better results are obtained when actual households are modelled, as we will see in part 2 of this tutorial.

The keys to defining this NetLogo model, so far, have been loading of a GIS map and various archaeologically established datasets, together number-crunching to calculate crop yields. When households are modelled, however, some code for human decision-making will also be needed. The model also illustrates how much functionality can be achieved with even a small NetLogo program!

See part 2 of this tutorial for the human side of the model. The full NetLogo program is at Modeling Commons.