GIS General Examples

extensions [ gis ]
globals [ cities-dataset
          rivers-dataset
          countries-dataset
          elevation-dataset ]
breed [ city-labels city-label ]
breed [ country-labels country-label ]
breed [ country-vertices country-vertex ]
breed [ river-labels river-label ]
patches-own [ population country-name elevation ]

to setup
  ca
  ; Note that setting the coordinate system here is optional, as
  ; long as all of your datasets use the same coordinate system.
  gis:load-coordinate-system (word "data/" projection ".prj")
  ; Load all of our datasets
  set cities-dataset gis:load-dataset "data/cities.shp"
  set rivers-dataset gis:load-dataset "data/rivers.shp"
  set countries-dataset gis:load-dataset "data/countries.shp"
  set elevation-dataset gis:load-dataset "data/world-elevation.asc"
  ; Set the world envelope to the union of all of our dataset's envelopes
  gis:set-world-envelope (gis:envelope-union-of (gis:envelope-of cities-dataset)
                                                (gis:envelope-of rivers-dataset)
                                                (gis:envelope-of countries-dataset)
                                                (gis:envelope-of elevation-dataset))
end 

; Drawing point data from a shapefile, and optionally loading the
; data into turtles, if label-cities is true

to display-cities
  ask city-labels [ die ]
  foreach gis:feature-list-of cities-dataset
  [ gis:set-drawing-color scale-color red (gis:property-value ? "POPULATION") 5000000 1000
    gis:fill ? 2.0
    if label-cities
    [ ; a feature in a point dataset may have multiple points, so we
      ; have a list of lists of points, which is why we need to use
      ; first twice here
      let location gis:location-of (first (first (gis:vertex-lists-of ?)))
      ; location will be an empty list if the point lies outside the
      ; bounds of the current NetLogo world, as defined by our current
      ; coordinate transformation
      if not empty? location
      [ create-city-labels 1
        [ set xcor item 0 location
          set ycor item 1 location
          set size 0
          set label gis:property-value ? "NAME" ] ] ] ]
end 

; Drawing polyline data from a shapefile, and optionally loading some
; of the data into turtles, if label-rivers is true

to display-rivers
  ask river-labels [ die ]
  gis:set-drawing-color blue
  gis:draw rivers-dataset 1
  if label-rivers
  [ foreach gis:feature-list-of rivers-dataset
    [ let centroid gis:location-of gis:centroid-of ?
      ; centroid will be an empty list if it lies outside the bounds
      ; of the current NetLogo world, as defined by our current GIS
      ; coordinate transformation
      if not empty? centroid
      [ create-river-labels 1
          [ set xcor item 0 centroid
            set ycor item 1 centroid
            set size 0
            set label gis:property-value ? "NAME" ] ] ] ]
end 

; Drawing polygon data from a shapefile, and optionally loading some
; of the data into turtles, if label-countries is true

to display-countries
  ask country-labels [ die ]
  gis:set-drawing-color white
  gis:draw countries-dataset 1
  if label-countries
  [ foreach gis:feature-list-of countries-dataset
    [ let centroid gis:location-of gis:centroid-of ?
      ; centroid will be an empty list if it lies outside the bounds
      ; of the current NetLogo world, as defined by our current GIS
      ; coordinate transformation
      if not empty? centroid
      [ create-country-labels 1
        [ set xcor item 0 centroid
          set ycor item 1 centroid
          set size 0
          set label gis:property-value ? "CNTRY_NAME" ] ] ] ]
end 

; Loading polygon data into turtles connected by links

to display-countries-using-links
  ask country-vertices [ die ]
  foreach gis:feature-list-of countries-dataset
  [ foreach gis:vertex-lists-of ?
    [ let previous-turtle nobody
      let first-turtle nobody
      ; By convention, the first and last coordinates of polygons
      ; in a shapefile are the same, so we don't create a turtle
      ; on the last vertex of the polygon
      foreach but-last ?
      [ let location gis:location-of ?
        ; location will be an empty list if it lies outside the
        ; bounds of the current NetLogo world, as defined by our
        ; current GIS coordinate transformation
        if not empty? location
        [ create-country-vertices 1
          [ set xcor item 0 location
            set ycor item 1 location
            ifelse previous-turtle = nobody
            [ set first-turtle self ]
            [ create-link-with previous-turtle ]
            set hidden? true
            set previous-turtle self ] ] ]
      ; Link the first turtle to the last turtle to close the polygon
      if first-turtle != nobody and first-turtle != previous-turtle
      [ ask first-turtle
        [ create-link-with previous-turtle ] ] ] ]
end 

; Using gis:intersecting to find the set of patches that intersects
; a given vector feature (in this case, a river).

to display-rivers-in-patches
  ask patches [ set pcolor black ]
  ask patches gis:intersecting rivers-dataset
  [ set pcolor cyan ]
end 

; Using gis:apply-coverage to copy values from a polygon dataset
; to a patch variable

to display-population-in-patches
  gis:apply-coverage countries-dataset "POP_CNTRY" population
  ask patches
  [ ifelse (population > 0)
    [ set pcolor scale-color red population 500000000 100000 ]
    [ set pcolor blue ] ]
end 

; Using find-one-of to find a particular VectorFeature, then using
; gis:intersects? to do something with all the features from another
; dataset that intersect that feature.

to draw-us-rivers-in-green
  let united-states gis:find-one-feature countries-dataset "CNTRY_NAME" "United States"
  gis:set-drawing-color green
  foreach gis:feature-list-of rivers-dataset
  [ if gis:intersects? ? united-states
    [ gis:draw ? 1 ] ]
end 

; Using find-greater-than to find a list of VectorFeatures by value.

to highlight-large-cities
  let united-states gis:find-one-feature countries-dataset "CNTRY_NAME" "United States"
  gis:set-drawing-color yellow
  foreach gis:find-greater-than cities-dataset "POPULATION" 10000000
  [ gis:draw ? 3 ]
end 

; Drawing a raster dataset to the NetLogo drawing layer, which sits
; on top of (and obscures) the patches.

to display-elevation
  gis:paint elevation-dataset 0
end 

to display-elevation-in-patches
  ; This is the preferred way of copying values from a raster dataset
  ; into a patch variable: in one step, using gis:apply-raster.
  gis:apply-raster elevation-dataset elevation
  ; Now, just to make sure it worked, we'll color each patch by its
  ; elevation value.
  let min-elevation gis:minimum-of elevation-dataset
  let max-elevation gis:maximum-of elevation-dataset
  ask patches
  [ ; note the use of the "<= 0 or >= 0" technique to filter out
    ; "not a number" values, as discussed in the documentation.
    if (elevation <= 0) or (elevation >= 0)
    [ set pcolor scale-color black elevation min-elevation max-elevation ] ]
end 

; This is a second way of copying values from a raster dataset into
; patches, by asking for a rectangular sample of the raster at each
; patch. This is somewhat slower, but it does produce smoother
; subsampling, which is desirable for some types of data.

to sample-elevation-with-patches
  let min-elevation gis:minimum-of elevation-dataset
  let max-elevation gis:maximum-of elevation-dataset
  ask patches
  [ set elevation gis:raster-sample elevation-dataset self
    if (elevation <= 0) or (elevation >= 0)
    [ set pcolor scale-color black elevation min-elevation max-elevation ] ]
end 

; This is an example of how to select a subset of a raster dataset
; whose size and shape matches the dimensions of the NetLogo world.
; It doesn't actually draw anything; it just modifies the coordinate
; transformation to line up patch boundaries with raster cell
; boundaries. You need to call one of the other commands after calling
; this one to see its effect.

to match-cells-to-patches
  gis:set-world-envelope gis:raster-world-envelope elevation-dataset 0 0
  cd
  ct
end 

; This command also demonstrates the technique of creating a new, empty
; raster dataset and filling it with values from a calculation.
;
; This command uses the gis:convolve primitive to compute the horizontal
; and vertical Sobel gradients of the elevation dataset, then combines
; them using the square root of the sum of their squares to compute an
; overall "image gradient". This is really more of an image-processing
; technique than a GIS technique, but I've included it here to show how
; it can be easily done using the GIS extension.

to display-gradient-in-patches
  let horizontal-gradient gis:convolve elevation-dataset 3 3 [ 1 0 -1 2 0 -2 1 0 -1 ] 1 1
  let vertical-gradient gis:convolve elevation-dataset 3 3 [ 1 2 1 0 0 0 -1 -2 -1 ] 1 1
  let gradient gis:create-raster gis:width-of elevation-dataset gis:height-of elevation-dataset gis:envelope-of elevation-dataset
  let x 0
  repeat (gis:width-of gradient)
  [ let y 0
    repeat (gis:height-of gradient)
    [ let gx gis:raster-value horizontal-gradient x y
      let gy gis:raster-value vertical-gradient x y
      if ((gx <= 0) or (gx >= 0)) and ((gy <= 0) or (gy >= 0))
      [ gis:set-raster-value gradient x y sqrt ((gx * gx) + (gy * gy)) ]
      set y y + 1 ]
    set x x + 1 ]
  let min-g gis:minimum-of gradient
  let max-g gis:maximum-of gradient
  gis:apply-raster gradient elevation
  ask patches
  [ if (elevation <= 0) or (elevation >= 0)
    [ set pcolor scale-color black elevation min-g max-g ] ]
end 


; Public Domain:
; To the extent possible under law, Uri Wilensky has waived all
; copyright and related or neighboring rights to this model.