Работа с растровыми геоданными в NetLogo

Импорт изображений

Для загрузки изображений в NetLogo и рисования линий с помощью черепашек существует особый слой игрового мира — слой для рисования (drawing layer), располагающийся между участками и движущимся по ним черепашками. Сразу после запуска NetLogo этот слой пуст и прозрачен.

Команда наблюдателя import-drawing позволяет импортировать файл с изображением с диска в слой для рисования. Поддерживаются следующие графические форматы: BMP, JPG, GIF и PNG.

import-drawing полезна для создания фона, изображения которое видно нам, но не заметно для черепашек или участков. Тем не менее, все, что мы нарисуем на экране при помощи черепашек, можно экспортировать с помощью меню File\Export\View или команды export-view “filename”, которая сохраняет текущее изображение игрового мира (всех его видимых слоев) во внешнем файле filename формата PNG (рекомендуется указывать имя файла filename с расширением .png).

Если нужно, чтобы черепашки или участки реагировали на изображение, используйте import-pcolors или import-pcolors-rgb.

import-pcolors “filename” cчитывает изображение из файла filename, изменяет его размер так, чтобы изображение поместилось на сетку из участков, сохраняя при этом пропорции изображения. Цвет пикселя изображения становится цветом соответствующего участка.

Однако цвета импортированного изображения могут быть изменены, т.к. цветовое пространство NetLogo не содержит всех возможных цветов. import-pcolors может работать медленно для больших изображений и при большом количестве участков.

Т.к. import-pcolors устанавливает значения the pcolor, то теперь агенты могут «чувствовать» изображение и взаимодействовать с ним.

Чтобы цвета импортированного изображения не искажались и агенты имели возможность взаимодействовать с импортированным изображением используется команда import-pcolors-rgb “filename”.

Для того чтобы изображения отображались в участки качественно, размер каждого участка нужно установить равным одному пикселю (Settings\Patch Size). Пример использования: Models Library\Code Examples\Image Import Examples

Очистка слоя для рисования

Чтобы очистить слой для рисования, используется команда clear-drawing, выполняемая наблюдателем. Вы также можете использовать clear-all, которая среди прочего, очищает и слой для рисования.

Рисование при помощи черепашек

Черепашки могут рисовать или удалять линии в игровом мире. Делают они это с помощью "пера", которое есть у каждой черепашки. Используя команду pen-down, черепашка может оставлять след "пера", соответствующий траектории ее движения. Команда pen-erase удаляет линии, лежащие на пути черепашки.

Работа с растровыми данными с помощью расширения GIS

Загрузка данных

Загрузка данных осуществляется аналогично векторным данным (с помощью gis:load-dataset, gis:envelope-of, gis:set-world-envelope и т.п.). Рассмотрим это на примере файла world-elevation.asc, где содержатся данные о высотах (elevation) объектов над уровнем моря (файл GIS Elevation0.nlogo из архива)

extensions [ gis ]
globals [ elevation-dataset ]

to setup
  clear-all
  set elevation-dataset gis:load-dataset "data/world-elevation.asc"
  gis:set-world-envelope (gis:envelope-of elevation-dataset)
end

; Рисуем растровые данные в слое для рисования, который находится над
; участками и незаметен для них.
to display-elevation
; Данные отображаются в оттенках серого (максимальное значение в наборе 
; данных - белый цвет, минимальное - черный). "0" определяет уровень 
; прозрачности (0 - непрозрачный, 255 - полностью прозрачный)
  gis:paint elevation-dataset 0
end

Сохранение растровых данных в переменных участков

Скопируем теперь данные о высотах из растрового набора и поместим в переменную участков elevation (GIS Elevation1.nlogo).

extensions [ gis ]
globals [ elevation-dataset ]
patches-own [ elevation ]

to setup
  clear-all
  set elevation-dataset gis:load-dataset "data/world-elevation.asc"
  gis:set-world-envelope (gis:envelope-of elevation-dataset)
end

; Рисуем растровые данные в слое для рисования.
to display-elevation
  gis:paint elevation-dataset 0
end

to display-elevation-in-patches
  ; Скопируем данные из растрового набора и сохраним значения высот (elevation) в
  ; переменной elevation.
  gis:apply-raster elevation-dataset elevation
  ; Окрасим каждый участок в цвет, соответствующий его высоте
  ; Определяем
  let min-elevation gis:minimum-of elevation-dataset ; минимальное значение в наборе
  let max-elevation gis:maximum-of elevation-dataset ; максимальное...
  ask patches
  [ ; Условие "<= 0 или >= 0" отфильтровывает значения "not a number"
    if (elevation <= 0) or (elevation >= 0)
    [ set pcolor scale-color patches-color elevation min-elevation max-elevation ] ]
end

Команда gis:apply-raster RasterDataset patch-variable копирует значения величины, хранящейся в растровом наборе данных RasterDataset, и сохраняет их в переменной, относящейся к участкам patch-variable. При необходимости передискретизирует (ресэмплирует) растр, так чтобы границы элементов растра соответствовали границам участков, составляющих игровой мир. Переменной участков, не покрытых растром, присваивается значение "not a number".

Команды gis:minimum-of RasterDataset, gis:maximum-of RasterDataset возвращают соответственно минимальное и максимальное значения в растровом наборе данных RasterDataset.

Альтернативный способ сохранения значений переменной, хранящейся в растре, в переменной участков (GIS Elevation2.nlogo)

to another-way-2-display-elevation
  let min-elevation gis:minimum-of elevation-dataset
  let max-elevation gis:maximum-of elevation-dataset
  ask patches
  [ ; raster-sample возвращает значение элемента растра (высоты) для указанного 
    ; расположения (каждого участка).
    set elevation gis:raster-sample elevation-dataset self
    if (elevation <= 0) or (elevation >= 0)
    [ set pcolor scale-color patches-color elevation min-elevation max-elevation ] ]
end

Команда gis:raster-sample RasterDataset sample-location возвращает значение величины, хранящейся в заданном месте (sample-location) растра. Расположение может быть:

• Списком длины 2: координаты точки в игровом мире ([ xcor ycor ]), вроде тех, что выдаются командой location-of Vertex. При необходимости, растр дискретизуется в точке с указанными координатами.
• Списком длины 4: «конверт» геоданных, возвращаемый командой envelope-of. При необходимости, растр дискретизуется на всей площади «конверта».
• Участком.
• Черепашкой.
• Вершиной (Vertex).

Если заданное расположение — точка, то дискретизация осуществляется с помощью метода, указанного в set-sampling-method. Если же расположение — площадь (конверт или участок), дискретизация осуществляется вычислением среднего значения для всех элементов растра, покрывающих заданную площадь.

Как и в gis:apply-raster, переменной участков, не покрытых растром, присваивается значение "not a number".

Замечание: gis:apply-raster работает быстрее, т.к. в нем осуществляется более простой метод передискретизации — метод ближайшего соседа. В качестве нового значения выбирается известное значение функции среди четырех соседних точек.

Выберем фрагмент растровых данных, размеры и форма которого совпадают с игровым миром (GIS Elevation3.nlogo).

to match-cells-to-patches
  gis:set-world-envelope gis:raster-world-envelope elevation-dataset x y
  cd
end

Основу здесь составляют две команды. * gis:raster-world-envelope RasterDataset x y — возвращает «конверт» для подмножества растрового набора данных RasterDataset, совпадающего размерами и формой с игровым миром (полученный «конверт» можно затем использовать в качестве аргумента для команд, осуществляющих преобразования, вроде set-transformation-ds). x и y — значения координат верхнего левого угла «конверта». Элементы растра нумеруются слева направо и сверху вниз, начиная с нуля. Поэтому верхний левый угол имеет координаты (0, 0), а нижний правый — (gis:width-of RasterDataset - 1, gis:height-of RasterDataset - 1). * gis:set-world-envelope gis-envelope, которая задает картографическое преобразование «конверта» геоданных gis-envelope в игровой мир NetLogo, при котором сохраняются одинаковые масштабы преобразований вдоль осей x и y.

Сама процедура match-cells-to-patches ничего не рисует. Для этого можно, например, воспользоваться одной из процедур рисования, приведенных выше.

Преобразование растрового набора данных

В качестве примера создадим пустой набор растровых данных и заполним его результатами расчетов: выделим из исходного набора контуры по методу Собеля (GIS Elevation4.nlogo).

; Создадим пустой набор растровых данных и заполним его результатами расчетов.
to display-gradient-in-patches
; Вычислим горизонтальный и вертикальный градиенты с помощью свертки (convolve)
; (выделение контуров по методу Собеля).
  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
  ; Создадим пустой набор данных gradient...
  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
      ; значениями общего "градиента изображения"
      gis:set-raster-value gradient x y sqrt ((gx * gx) + (gy * gy))
      ; переходим к новым элементам из gradient
      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

gis:convolve RasterDataset kernel-rows kernel-columns kernel key-column key-row — возвращает новый растр, данные которого состоят из данных заданного растра RasterDataset, свернутого (convolved) с заданным ядром (kernel).

Свертка в компьютерной графике — это операция вычисления нового значения выбранного пикселя (в нашем случае — элемента растра), учитывающая значения окружающих его пикселей. Для вычисления значения используется матрица, называемая ядром свертки. Обычно ядро свертки является квадратной матрицей n х n, где n — нечетное, однако ничто не мешает сделать матрицу прямоугольной. Во время вычисления нового значения выбранного пикселя ядро свертки как бы «прикладывается» своим центром (именно тут важна нечетность размера матрицы) к данному пикселю. Окружающие пиксели так же накрываются ядром. Далее высчитывается сумма, где слагаемыми являются произведения значений пикселей на значения ячейки ядра, накрывшей данный пиксель. Сумма делится на сумму всех элементов ядра свертки. Полученное значение как раз и является новым значением выбранного пикселя. Если применить свертку к каждому пикселю изображения, то в результате получится некий эффект, зависящий от выбранного ядра свертки.

• kernel-rows, kernel-columns — число строк и столбцов ядра свертки;
• key-column, key-row — номер столбца и строки ключевого (центрального) элемента ядра;
• kernel — матрица-ядро свертки.

Значения элементов ядра представляются виде списка, в котором элементы матрицы представлены слева направо и сверху вниз. Для матриц 3х3 это выглядит следующим образом:

Например, для ядра оператора Собеля:

можно записать следующее (через key обозначен ключевой элемент):

let horizontal-gradient gis:convolve dataset 3 3 [ 1 0 -1 2 0 -2 1 0 -1 ] 1 1

Здесь

• 3 3 — размерность матрицы-ядра;
• [ 1 0 -1 2 0 -2 1 0 -1 ] — ядро;
• 1 1 — расположение ключа.

Напомним, что нумерация в NetLogo начинается с 0.

Результат применения метод Собеля (display-gradient-in-patches) к тестовому изображению (вверху) приведен ниже.

Взаимодействие агентов с растром

В GIS Downhill.nlogo моделируется как «капли» дождя стекают по склонам каньона и скапливаются во впадинах.

extensions [ gis ]
globals [ elevation ]
patches-own [ friction ]

to setup
  clear-all
  set elevation gis:load-dataset "data/local-elevation.asc"
  gis:set-world-envelope gis:envelope-of elevation
  gis:apply-raster elevation friction
  ask patches
  [ sprout 1 [ set color blue ]]
  gis:paint elevation 0
end

to go
  ask turtles [ downhill friction ]
  tick
end

Ключевое значение здесь имеет команда downhill, формат которой: downhill patch-variable. Она перемещает черепашку в один из соседних участков (всего их 8) с наименьшим значением переменной участков patch-variable. Если значения переменной в соседних участках больше, чем в текущем, то черепашка остается на месте. Если среди соседних участков есть несколько с наименьшими значениями, черепашка выбирает один из них случайным образом. Нечисловые значения переменной patch-variable игнорируются.

sprout число [ команды ] — создает заданное число новых черепашек на текущем участке и приказывает им выполнить команды из блока.

На рис. ниже показано начало (вверху) движения «капель» и положение после нескольких десятков шагов (внизу).

Более продвинутая модель движения капель показана ниже (GIS Gradient.nlogo).

extensions [ gis ]
globals [ elevation slope aspect ]

to setup
  clear-all
  set elevation gis:load-dataset "data/local-elevation.asc"
  gis:set-world-envelope gis:envelope-of elevation
  ; вычисляем градиенты (наклоны) в вертикальном и горизонтальном направлении
  let horizontal-gradient gis:convolve elevation 3 3 [ 1 0 -1 1 0 -1 1 0 -1 ] 1 1
  let vertical-gradient gis:convolve elevation 3 3 [ 1 1 1 0 0 0 -1 -1 -1 ] 1 1
  ; создаем два растра: slope – для хранения величины градиента в данной точке и
  ; aspect – для ориентации черепашки
  set slope gis:create-raster gis:width-of elevation gis:height-of elevation gis:envelope-of elevation
  set aspect gis:create-raster gis:width-of elevation gis:height-of elevation gis:envelope-of elevation
  let x 0
  repeat (gis:width-of slope)
  [ let y 0
    repeat (gis:height-of slope)
    [ ; узнаем величину градиентов в данной точке
      let gx gis:raster-value horizontal-gradient x y
      let gy gis:raster-value vertical-gradient x y
      ; вычисляем полный градиент
      gis:set-raster-value slope x y sqrt ((gx * gx) + (gy * gy))
      ifelse (gx != 0) or (gy != 0)
      ; черепашка должна ориентироваться по направлению градиента
      [ gis:set-raster-value aspect x y atan gx gy ]
      [ gis:set-raster-value aspect x y 0 ]
      set y y + 1 ]
    set x x + 1 ]
  ; для отображения растра используется билинейная интерполяция
  gis:set-sampling-method aspect "bilinear"
  ask patches ; для всех участков
  [ sprout 1 ; создадим в каждом по черепашке
    [ set color blue ; синего цвета
      let h gis:raster-sample aspect self ; h хранит ориентацию черепашки
      ifelse h >= -360
      [ set heading subtract-headings h 180 ] 
      [ die ] ] ]
  gis:paint elevation 0
end

to go
  ask turtles
  [ ; Каждая черепашка перемещается на случайное расстояние вдоль направления aspect
    forward random-normal 0.1 0.1
    let h gis:raster-sample aspect self
    ifelse h >= -360
    [ set heading subtract-headings h 180 ]
    [ die ] ]
  tick
  if not any? turtles ; закончить работу, если не осталось черепашек
  [ stop ]
end

Несколько команд требуют дополнительных пояснений.

• atan x y превращает смещения x и y в ориентацию черепашки, в градусах (от 0 до 360). Заметим, что реализация atan рассчитана на геометрию игрового мира NetLogo, где направление 0 означает ориентацию черепашки вертикально вверх, 90 — вправо, и так далее по часовой стрелке (обычно в геометрии угол 0 означает направление вправо, 90 — вверх, а движение идет против часовой стрелки). Когда y = 0: если x > 0, atan возвращает 90; если x < 0, то atan возвращает 270; если x = 0 atan выдает сообщение об ошибке.
• subtract-headings heading1 heading2 вычисляет разность между двумя заданными ориентациями (headings), т.е. значение в градусах наименьшего угла на который нужно повернуть heading2, чтобы совместить его с heading1. Положительный результат означает вращение по часовой стрелке, отрицательный — против часовой. Результат всегда находится в интервале (-180; 180].

show subtract-headings 80 60
=> 20
show subtract-headings 60 80
=> -20
show subtract-headings 5 355
=> 10
show subtract-headings 355 5
=> -10
show subtract-headings 180 0
=> 180
show subtract-headings 0 180
=> 180

• heading — встроенная переменная черепашек, указывающая в каком направлении ориентирована черепашка. Значения находятся в интервале [0; 360). 0 — север, 90 — восток, и т. д.

Примеры:

set heading 45      ;; черепашка ориентирована на северо-восток
set heading heading + 10 ;; что же, что "rt 10"

На рис. показаны результаты работы модели после 23 (вверху) и 127 (внизу) шагов.

В следующей программе (GIS Canyon.nlogo) сделано сразу несколько улучшений:

• Использован новый род агентов — «капли» (raindrops).
• Улучшено отображение цветов.
• Пользователь может добавлять «капли».
• «Капли» могут оставлять следы на поверхности.
• Контролируемая интенсивность осадков (rain-rate).
• Направление движения капли зависит от высоты поверхности и толщины слоя воды на ней.

Вот что получается в результате:

А вот как это реализовано:

extensions [ gis ]
breed [ raindrops raindrop ]

patches-own [ elevation ]

globals [
  elevation-dataset
  color-min
  color-max
  water-height    ;; высота одной единицы воды
  border          ;; множество граничных участков
]

to startup
; считывание геоданных удобно организовать еще при открытии модели
  set elevation-dataset gis:load-dataset "data/Grand Canyon data.asc"     
  gis:set-world-envelope (gis:envelope-of elevation-dataset)   
  gis:apply-raster elevation-dataset elevation
  ; добавим немного высоты, чтобы самые высокие точки не отображались чисто белыми
  set color-max gis:maximum-of elevation-dataset + 200
  let min-elevation gis:minimum-of elevation-dataset
  ; минимумы также не должны быть черными
  set color-min min-elevation - ((color-max - min-elevation) / 10) 
  ; копируем растровые данные в участки
  set-default-shape turtles "circle"
  setup
end

;; удалим "капли" и следы их течения, оставшиеся от прошлых запусков
;; и зададим некоторые глобальные переменные
to setup
  clear-drawing
  clear-turtles
  ask patches
  [ if (elevation <= 0) or (elevation >= 0) ; проверка на "not a number"
    [ set pcolor scale-color brown elevation color-min color-max ] ] 
      ; отображаем участки 
  set water-height 10
  set border patches with [ count neighbors != 8 ] ; граничные участки те, что не имеют достаточно соседей
end

to go
  ; Создаем новые "капли" кликаньем левой кнопкой мыши.
  ; Проверяем, пришелся ли клик на пустой участок
  if mouse-down? and not any? turtles-on patch mouse-xcor mouse-ycor
  [ ; Если да, то создаем на этом участки "каплю" красного цвета 
    create-raindrops 1 ; создаем черепашку, являющуюся "каплей"...
    [ setxy mouse-xcor mouse-ycor ; там, где кликнули
      set size 2
      set color red
    ]
  ]
  ; остальные "капли" падают случайным образом
  ; интенсивность "дождя" задается rain-rate
  create-raindrops rain-rate
  [ move-to one-of patches
    set size 2
    set color blue ]

  ifelse draw?
    [ ask turtles [ pen-down ] ]
    [ ask turtles [ pen-up ] ]

  ask raindrops [ flow ]

  ask border
  [ ; "капли", утекающие за границу мира, уничтожаются
    ask turtles-here [ die ]  ]
  tick
end

to flow ; выполняют черепашки
  ; Высота участка складывается из высоты поверхности и толщины слоя воды на ней.
  ; Из соседних участков выбирается участок с наименьшей высотой - target
  let target min-one-of neighbors [ elevation + (count turtles-here * water-height) ]
  ; если текущего участка имеет высоту больше чем target, то перемещаемся на target
  if (elevation + (count turtles-here * water-height)) >
     [elevation + (count turtles-here * water-height)] of target
    [ move-to target ]
end

Можно добавить в программу простейшую модель эрозии (GIS Canyon Eroded.nlogo). Для этого вводятся глобальные переменные erosion-rate (темп эрозии) и soil-amount (количество грунта, уносимого одной черепашкой):

to erode
  if random-float 1 < erosion-rate ; если есть возможность для эрозии
  [ let washout soil-amount * count turtles-here 
    ; то количество уносимого с участка грунта
    set elevation elevation - washout            
    ; пропорционально числу "капель" на нем 
    if (elevation <= 0) or (elevation >= 0)
    [ set pcolor scale-color brown elevation color-min color-max ]
  ]
end

Однако здесь можно столкнуться с проблемой: «капля» сама вырывает себе яму, из которой уже не может выбраться. Необходимо подбирать правильное соотношение между soil-amount и water-height. К тому же работа программы существенно замедляется из-за перерисовки рельефа.

Преобразование изображений в формат .asc и .grd

Растровые данные можно получить, в частности, конвертируя рисунки распространенных графических форматов (.bmp, .jpg, .gif, .png и других) с помощью свободно распространяемой ГИС LandSerf.

Коллекции растровых данных в Internet

Данные о рельефе США можно найти в виде National Elevation Data (NED) на сайте Геологической службы США (USGS Seamless Data Distribution website) по адресу http://seamless.usgs.gov/. Скачанный файл приводится в порядок и сохраняется в .asc или .grd с помощью LandSerf.

Разное

Как вычислить минимальные и максимальные координаты участков?

show min-pxcor
show max-pxcor
show min-pycor
show max-pycor

Протяженность игрового мира «в участках» устанавливается кнопкой "Settings"

Как определить количество столбцов и строк в растре?

show gis:width-of rasterdataset
show gis:height-of rasterdataset

Как отметить расположение объекта с помощью цвета?

Допустим, интересующий нас объект находится в участке с координатами [20 30]. После того, как набор растровых данных спроектирован в игровой мир, набираем:

clear-all
ask patch 20 30 [set pcolor red]

Как узнать координаты участка?

Навести на него указатель мыши, и кликнуть правой кнопкой. В появившемся меню будут указаны интересующие координаты.

Как создать черепашку заданного цвета в заданном месте?

Создаем красную черепашку на участке с координатами [50 50]

ask patch 50 50 [sprout 1 [set color red]]

Как начертить траекторию движения черепашки?

Опустить перо:

ask turtles [pen-down]

Как передать значение переменной, расположенной в заданном месте растра в переменную участка, расположенную в том же месте?

patches-own [friction]
ask patch 50 50 [set friction gis:raster-value elevation 50 50]

Проверим:

show [friction] of patch 50 50
show gis:raster-value elevation 50 50

Как скопировать значение переменной из растрового набора данных в переменную, относящуюся к участкам?

gis:apply-raster dataset patch-variable

Количество участков в направлениях осей x и y и количество строк и столбцов растра могут не совпадать. При этом автоматически осуществляется передискретизация (resampling) растра. Для того, чтобы качество изображения в игровом мире отвечало качеству растровых данных, размер сетки из участков должен совпадать с размерами растра (количеством строк и столбцов).