-
Notifications
You must be signed in to change notification settings - Fork 13
MapAlgebra
This document is under development.
Pedro R. Andrade
It is possible to implement map algebra procedures using some TerraME concepts. This tutorial describes the basic structure of a script that executes different types of map algebra operations.
Map Algebra (From Wikipedia)
Map algebra is a set-based algebra for manipulating geographic data, proposed by Dr. Dana Tomlin in the early 1980s. It is a set of primitive operations in a geographic information system (GIS) which allows two or more raster layers ("maps") of similar dimensions to produce a new raster layer (map) using algebraic operations such as addition, subtraction etc.
Depending on the spatial neighborhood, GIS transformations are categorized into four classes: local, focal, global, and zonal. Local operations works on individual raster cells, or pixels. Focal operations work on cells and their neighbors, whereas global operations work on the entire layer. Finally, zonal operations work on areas of cells that share the same value. The input and output for each operator being map, the operators can be combined into a procedure or script, to perform complex tasks.
To implement a script that uses map algebra operations, it is necessary to understand the types and functions below for the different types of operations. If you are not familiar with these concepts, it is recommended to read An Introduction to TerraME and Creating and Filling Cellular Spaces in TerraME.
Type or function | Local | Focal | Global | Zonal |
---|---|---|---|---|
Project |
X | X | X | X |
CellularSpace |
X | X | X | X |
forEachCell() |
X | X | X | X |
forEachCellPair() |
X | X | X | X |
save() |
X | X | X | X |
Map |
X | X | X | X |
createNeighborhood() |
X | X | X | |
forEachNeighborhood() |
X | X | X | |
split() |
X | |||
Trajectory |
X |
Local operations use existing attribute values to compute new attributes or to update existing ones.
To work with Map Algebra in TerraME, the data needs to belong to a Project
. If it already belongs
to a Project
, this step can be skipped.
import("gis")
proj = Project{
file = "amazonia.tview",
clean = true,
amazonia = filePath("amazonia.shp")
}
Once the data already belongs to a Project
, it is necessary to load the data as a CellularSpace
.
cs = CellularSpace{
project = proj,
layer = "amazonia"
}
After that, it is now possible to execute map algebra operations using forEachCell
, which
takes a function as argument and applies it over every Cell
of the CellularSpace
. The code
below creates attribute distweight
by inverting an existent attribute called distroad
.
At this point, it is possible to create as many attributes as needed using basic Lua
facilities.
forEachCell(cs, function(cell)
cell.distweight = 1 / cell.distroad
end)
The data created can be saved into another shapefile using save()
. The first argument
gets the new Layer
's name (that will also be the file name to be saved), and the
second argument gets the new attribute(s) to be saved. This function also saves all
the original attributes in the file.
cs:save("myamazonia", "distweight")
Now it is possible to load the CellularSpace
using only argument file
.
Alternatively, it is possible to use project
and layer
arguments instead of file
.
cs = CellularSpace{
file = filePath("myamazonia.shp")
}
Map{
target = cs,
select = "distweight",
min = 0,
max = 1,
slices = 8,
color = "Reds",
invert = true
}
An example to rename attributes is shown below.
The file amazonia.shp
from base package was created from the previous version of TerraLib (4)
and now needs to be updated to the current version (5). First we need to create a Project
with a Layer
containing the file. See the code below.
import("gis")
proj = Project{
file = "amazonia.tview",
clean = true,
amazonia = filePath("amazonia.shp")
}
In the previous version of TerraLib, a cellular space had the attributes Lin
and Col
.
In version 5, the attributes are now
row
and col
. Additionally, in the previous version the fist line of the cellular space started from the bottom
and now it starts from the top. To consider these differences, we need to set the arguments xy
and zero
while loading a CellularSpace
.
The attributes x
and y
of the CellularSpace
are then created using this configuration.
cs = CellularSpace{
project = proj,
layer = "amazonia",
zero = "bottom",
xy = {"Col", "Lin"},
}
After loading the data, we can now create row
and col
from y
and x
, respectively.
As we want to rename an attribute from Col
to col
, it is necessary to set the
first one as nil
because they have the same lower case, and TerraLib requires that
attributes in a cellular space must have different lower cases.
forEachCell(cs, function(cell)
cell.row = cell.y
cell.col = cell.x
cell.Col = nil
end)
After that, we only need to save()
the CellularSpace
into a file.
The two new attributes need to be passed as second argument to save.
cs:save("myamazonia", {"row", "col"})
Now it is possible to load the CellularSpace
using only argument file
.
Neighborhoods and Maps now work properly.
cs = CellularSpace{
file = filePath("myamazonia.shp")
}
cs:createNeighborhood()
Map{
target = cs,
select = "defor_10",
min = 0,
max = 1,
slices = 8,
color = "RdYlGn",
invert = true
}
If you have comments, doubts or suggestions related to this document, please write a feedback to pedro.andrade <at> inpe.br.
Back to wiki or terrame.org.