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:
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
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
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
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
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
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:
Pingback: NetLogo post summary | Scientific Gems
Pingback: Revisiting Artificial Anasazi: a NetLogo tutorial (part 1) | Scientific Gems
Thank you for your information. i am trying to make evacuation model by using Letlogo. Before seeing your article, I just used ‘import-pcolor’. And I realize that it must need GIS data to make better agent-based model.
Is this an evacuation model for a building? You wouldn’t really need GIS data for that.
I am working on an evacuation model for a building. I used shape file and vectordataset for setting up the environment. What is your suggestion for generating the environment?
Also, I got a problem with the gis:apply coverage primitive. It gives a runtime error. “Extension exception: not a vectordataset ”
Do you know any solution for this?
I can send you both shap file and the netlogo code, if you need.
The error “Extension exception: not a vectordataset” may mean that you have a raster file rather than a vector file. As to generating the environment, it depends on your goals for the simulation.
May I know, how did the avg patch area of 510,000,00/N square KM calculated?
510,000,000 sq km is the surface area of the planet, and N is the number of patches. Of course, this has substantial distortion due to the map projection, as noted in the blog post.
I have lat/long data to map the lower 48 states of the U.S. Plus lat/long for selected cities. However, I wanted to add Congressional Districts or maybe zip code areas to the map, so I am trying the GIS extension to Netlogo. The maps look much better. However, when I try to put the cities on the vector map, they are off. Not sure how to coordinate the two. The Netlogo world sets 0,0 in the lower left corner, and I convert the city lat/long to Netlogo coordinates. Also do this for the non-GIS version of the maps, and that looks okay. But when I use the GIS map, it’s not fitting into the Netlogo world so that the lat/long data matches. Any way to fix this? These cities are not going to be in a GIS database – or at least many of them won’t. Thanks for any advice you have.
Well, first, there’s an option when you right-click edit on the Netlogo world to set the origin 0,0 (I usually have it in the centre).
Second, the Netlogo coordinates specify patch centres. If you have the origin in the lower left, x actually runs from -0.5 to max-pxcor + 0.5, and y runs from -0.5 to max-pycor + 0.5 (although only the range 0 to max-pxcor and 0 to max-pycor are legal turtle positions). You need to do the conversion correctly.
Third, there may still be a bug in the Netlogo GIS extension. I notice that the patches and lines in my example don’t 100% line up. You might need to do some experimentation with the coordinates of geographic features and then compensate for the error. I must confess that I’ve never done anything requiring that level of accuracy.
I think your second comment may help – trial-and-error resulted in “correction factors” very close to the min latitude and longitude that is subtraced from each city’s location. These are within 0.5 of the actual boundary lat/long points. Will have to look at this when I’m more awake.
The CONUS map looks pretty good, though. Cities seem to be where they should be – i.e., not yet under water and in the expected states! Need to check this out for individual states, though.
Glad to hear that you seem to have gotten a handle on what is happening.
Some more progress, but more questions. I imported the data files I use in Netlogo into QGIS. The cities lined up really nicely. Then I created a shapefile and went back to Netlogo. I can now create a U.S. map with the cities and have a button that turns the cities on and off. So then I added the cities as agents using the same data file that QGIS used. And they are all a little off. And it’s not consistently to one direction – some to the right, some up to the right, some down, etc. So is there a way to create or link an agent to a point from the GIS city? I am just getting into the GIS extension, so I am not familiar with all the options. Or is it possible to retrieve the patch from the GIS city points? That might work, too. Thanks for your help.
It’s hard to know what’s wrong without seeing your code. It could still be the kind of coordinate mismatch I described before (especially if the errors go the same way in similar parts of the map). Whatever it is, you could probably reduce the problem by editing the world (right-click and edit) to have more and smaller patches.
I don’t understand your other questions, precisely, You can certainly use the ordinary patch-here to identify the patch that an agent is on and setxy to reset an agent’s position.
Turns out that I can import the cities file into QGIS and create a shapefile. Then back in Netlogo the cities and background maps can be merged. The code example in the Models library that you mention in the tutorial shows how to do that and create agents. So all is well! Thanks for your help.
I’m trying to include 5 rasters in my model, all the rasters are from the same region, but when I apply more than 1 occurs the following extension error.
Extension exception: Dimensions (width=334991 height=334992) are too large
error while observer running GIS:APPLY-RASTER
called by procedure DISPLAY-CLASES-IN-PATCHES
called by Button ‘display-clases-in-patches’
What can I do? I think is something related with
gis:set-world-envelope (gis:envelope-union-of (gis:envelope-of apt-dataset) (gis:envelope-of clases-dataset))
Is something wrong?
You should probably contact the NetLogo people themselves, but I’m guessing your rasters are just much to large (if the dimensions are really width=334991 height=334992). Alternatively, if the rasters are from different areas, NetLogo is being forced to build such a large raster internally to hold the five pieces.
Can you combine them all the rasters into one using GIS software before loading them into NetLogo? It sounds like you may need to downsample the rasters to fit your patch grid as well.
Thank you very much, your info was very useful.
The problem was in the projections of the rasters and, as you said, NetLogo was being forced to build a large raster.
I loaded a shapefile of rivers into Netlogo and I want to determine for every patch the distance to the closest river. I thought I should first tell a patch whether it belongs to a river, but I can’t quite figure that out using gis:apply-coverage. One reason is that I don’t know how to access the attributes from a shapefile in Netlogo. Do you have some suggestions on how to do this?
The latest version of the GIS extension (see here) has some operators that look like they will do the job. Unfortunately patches can’t interact with a drawing of the rivers. You could, however, use gis:apply-coverage (which copies attributes from a shapefile) if there are attributes properties, or you can use gis:intersects?.