11. 栅格数据处理¶
[rastfunk]
从简单的查询到复杂的代数功能和逻辑运算，实际上，任何分析和建模在GRASS中都能够完成。
由于现有的数据数据分析模块太多(超过100个)，所以本章只对栅格数据处理做一个概括的介绍。关于遥感，涉及到卫星和航空影像的处理和分析，不是本书的内容。 涉及遥感分析的文献在参考书目已经列出。
正如第[sec:grassgui]章中描述的那样，大多数的分析可以用图形用户界面TclTkGRASS和GIS管理器来执行。这里我们将对每一个分析步骤做GUI操作和命令行操作的解释。为了熟悉每个模块、它们的用法和它们所需的参数，我们强烈推荐您使用命令行方式。
栅格值以它的位置坐标x,y(像元中心)和z值来定义，z值对应一个测量值或对象值(通常赋一个显示时使用的颜色值或灰度值)。有两种基本的操作可以用来生成栅格数据的专题图：
- 逐点的、面向栅格像元或像素的操作 (邻近操作)
- 面向矩阵或像素窗口的操作 (移动窗口操作)
除了”现有”的GRASS模块(提供一个特定的操作)，我们还能定义面向解决方案的操作，并通过算术模块来应用(请参阅第[sec:mapcalc]章)。
对于栅格数据各组成部分(像地图的空间参考、属性和像素的颜色值)的基本管理都由各自的栅格模块直接操作，因此，以前版本(GRASS 5.0和5.3)的模块不再需要。此外，创建地图统计值的模块，也只对选中地图的显示区域有效。
几乎所有的GRASS模块都有帮助文件，里面描述了模块功能并解释了命令的语法。在终端输入模块命令并使用参数，可以显示简要的帮助：
用 module name可以显示详细的帮助，里面有模块的描述和示例，这与GRASS主页上的帮助页面是一样的。使用该命令时，会自动打开一个浏览器来显示相应的帮助文件：
11.1. 显示栅格地图¶
[visu]
栅格数据的显示通过模块来完成，地图会显示在所谓的X监视器中。GRASS可以平行操作7个不同的监视器(x0,x1,x2,..,x6)：
d.mon x0
d.rast
用可以得到某一范围更为详尽的视图：
用可以定义某一特殊裁减区域或整个工作区域。分辨率为10米的默认区域可以这样设置：
g.region -d res=10.0 -pa
d.erase
d.rast rastermap
注意，用更改的设置必须通过传递给监视器才会起作用。也可以直接从栅格地图来调整范围和分辨率：
g.region rast=rastermap -p
d.redraw
d.rast rastermap1
d.rast -o rastermap2
只有当叠置地图rastermap2的像素值为NULL值(表示数据的间隙)时，下层底图rastermap1的相应像素才可见。也可以叠置地图 – 如对Spearfish示例数据中TK24的查验。
用下面的命令可以显示栅格地图相应的图例：
11.2. 查询栅格像元值和元数据¶
[rastquest]
在显示了栅格地图以后，为了得到像素的坐标和栅格属性，可以使用：
它也可以查询多幅没有显示的栅格地图：
用鼠标点击一个像素，它的特征值会显示在终端里。也能够查询多幅显示地图的同一像素的属性值。点击右键退出该模块。
用来显示栅格地图的基本信息和元数据。它同时会显示数据描述、数据类型、地图投影和类别的值域。
r.info landuse
r.info -r landuse
用可以控制栅格地图的属性，创建一个带标签数值和相应文本属性的表格：
r.cats map=landuse
1 residential
2 commercial and services
3 industrial
4 other urban
5 reservoirs
6 bare exposed rock
7 quarries, strip mines and gravel pits
8 transportation and utilities
为了确定Spearfish中各个地质单元的范围，可以使用如下显示的(请参阅)。该模块根据当前范围和栅格分辨率计算统计值，因此，如果在运行该命令前进行了缩放，那么将只计算缩放后的区域：
g.region rast=geology -p
r.report -h geology units=h
+---------------------------------------------------------------------+
| Category Information | |
|#|description | hectares|
|---------------------------------------------------------------------|
|1|metamorphic. . . . . . . . . . . . . . . . . . . . . . .| 1051.000|
|2|transition . . . . . . . . . . . . . . . . . . . . . . .| 13.000|
|3|igneous. . . . . . . . . . . . . . . . . . . . . . . . .| 3285.000|
|4|sandstone. . . . . . . . . . . . . . . . . . . . . . . .| 6755.000|
|5|limestone. . . . . . . . . . . . . . . . . . . . . . . .| 5537.000|
|6|shale. . . . . . . . . . . . . . . . . . . . . . . . . .| 4170.000|
|7|sandy shale. . . . . . . . . . . . . . . . . . . . . . .| 1019.000|
|8|claysand . . . . . . . . . . . . . . . . . . . . . . . .| 1307.000|
|9|sand . . . . . . . . . . . . . . . . . . . . . . . . . .| 3295.000|
|*|no data. . . . . . . . . . . . . . . . . . . . . . . . .| 168.000|
|---------------------------------------------------------------------|
|TOTAL |26,600.000|
+---------------------------------------------------------------------+
要处理一个时间序列的地图，通常需要给地图一个时间戳(地图创建时间，数据获取时间等)。可以完成这一任务。可以使用绝对和相对日期，它们独立地存储在文件历史信息中：
r.timestamp landuse date="27 Sep 2003"
r.timestamp landuse date="27 Sep 2003/20 Feb 2004"
用可以找到关于时间戳格式更为详细的描述。
11.3. 不同的栅格应用¶
[rastanwend]
GRASS对于栅格的分析有十分广泛的应用。本节我们将介绍一些非常常用的GRASS栅格模块。这里讨论的模块和应用只是很少的一部分。
11.3.1. 横断面分析¶
[transekt]
为了在栅格图表上显示一个横断面，可以使用菜单驱动的模块(见图[fig:dprofile])。
用能够将沿一个或多个横断面的、栅格像元或中值及算数平均值以ASCII格式显示出来。横断面能够以坐标对给出，也可以交互式地选择起点和终点(使用选项)。为提供了一个前端界面。它只需要横截面的起点，终点会根据和参数来计算。
r.profile
r.transect
[tb!]
[fig:dprofile]
11.3.2. 视线分析¶
[rlos]
可以基于高程地图执行视线分析。起点、此点距地面的高程以及视线分析的距离都可以通过坐标指明。
作为示例，我们将使用Spearfish数据集中的高程地图。本例中的初始坐标已经指定，但也可以用来确定：
# Put the setting on relief model but 20m resolution
g.region rast=elevation.10m res=20 -pa
# Calculation of visibility of a tower 15m over terrain
# (computionally extensive at 10m original resolution...)
r.los in=elevation.10m out=visibility coord=593670,4926877 \
obs=15 max=30000
# Control
d.erase
r.shaded.relief elevation.10m units=meters
d.rast elevation.10m_shade
d.rast -o roads
d.rast -o visibility
11.3.3. 地图叠置¶
[rmergepatch]
正如前面提到的那样，我们能够将栅格地图叠置显示，然而模块能够将叠置地图存储为一幅新的地图。所以，通过该模块我们可以将多幅地图合并为一幅地图。
第一幅地图在顶上，其他地图相继出现在下面，并且只有当上层地图的值为时，下层的地图才会显示。例如，Spearfish数据集中道路和地质图的叠置。需要注意的是，GRASS模块基本上工作在当前的范围和分辨率下。因此，区域的分辨率和范围必须提前调整好。
g.region rast=geology,roads -p [res=12.5] [-a]
r.patch in=roads,geology out=roads.on.geol
如果地图的顺序改变了，那么整个道路图层会被地质图覆盖。这样的输出地图基本上和原来的地质图是一样的。
11.3.4. 栅格数据的缓冲区分析¶
[rbuffer]
允许用户对栅格数据做缓冲区。该功能可以为每个道路类别创建噪音保护区，使用Spearfish数据集，首先列出已有的类别：
# Which road categories are available?
g.region rast=roads -p
r.report -h roads
|1|interstate
|2|primary highway, hard surface
|3|secondary highway, hard surface
|4|light-duty road, improved surface
|5|unimproved road
|*|no data
问题：只为州际公路(interstate)定义距离为100、250和500米的缓冲区。初始地图的分辨率是30米，练习中始终要注意这一点。
# Extracting the category 'interstate' with r.reclass
r.reclass roads out=interstate << EOF
1 = 1 interstate
EOF
# or extracting the category 'interstate' with r.mapcalc
r.mapcalc "interstate=if(roads==1,roads,zero())"
# Creating buffer zones
r.buffer in=interstate out=interstate.buf dist=100,250,500
# Control
d.rast.leg interstate.buf
[H]
[abb:strpuff]
得到的地图显示了有三个预定义缓冲区的州际公路。指定的缓冲区宽度总是从被缓冲要素(本例是州际公路)的中间起算，并且以区域定义的地图单位(如：米)来计算。地图单位可以用来查询：
g.proj -p
...
-PROJ_UNITS------------------------------------------------
unit : meter
units : meters
meters : 1.0
栅格地图转换为矢量数据可以在 [rtovect]中找到。
11.4. 色彩映射表的修改和分配¶
[farbtabellen]
影像文件中的值在GRASS中是通过色彩映射表来编码的。创建一个栅格影像后，栅格图便被分配一个默认的色彩映射表“rainbow”，很少有例外。有几种途径可以来创建特定地图的色彩映射表。
用可以分配一种标准的色彩映射表：
r.colors map=raster map color=standard table
r.colors map=raster map color=special table
color options: aspect,grey,grey.eq,grey.log,byg,byr,gyr,rainbow,ramp,
random,ryg,wave,rules
rules options: aspect,bcyr,byg,byr,elevation,evi,grey,gyr,rainbow,ramp,
ryg,slope,srtm,terrain,wave
指明变量，我们可以分配自己的色彩映射表：
r.colors map=geology color=rules << EOF
4 100 200 0
5 255 130 7
6 100 129 187
7 222 180 39
9 43 18 200
EOF
为了将色彩映射表从一个栅格图转移到另一个栅格图，或者是分配一个单独的或已经指定的色彩映射表，我们使用模块及参数”rast”：
最新分配的色彩映射表可以用来显示：
11.5. 地图的统计值¶
[statistik]
GRASS提供了内建的模块来计算地图的统计值。此外，对于更为复杂的地统计分析(见参考书目)，有用于统计软件R的接口。
直方图和统计值对于影像的增强和分析(对比度的调整、拉伸)是十分重要的信息。影像的直方图可以用来生成：
在交互式模式下，像素的分布表现为条(bar)，也可以显示为饼图。
模块将像素的分布表现为表格。该模块能被其他程序扩展，因此在使用前请先查看帮助页面。例如，可以同时查询多幅栅格地图：
在这里我们可以方便地使用模块，它是的封装。该模块能够交互地控制像素的分布、范围的分析等，并显示为表格。例如土地利用和地质的查询：
g.region rast=landuse -p
r.report -hen landuse,geology units=h
+---------------------------------------------------------------------+
| Category Information | |
|#|description | hectares|
|---------------------------------------------------------------------|
|1|residential | 676.00000|
| |--------------------------------------------------------|----------|
| |1|metamorphic. . . . . . . . . . . . . . . . . . . . . .| 23.00000|
| |3|igneous. . . . . . . . . . . . . . . . . . . . . . . .| 18.00000|
| |4|sandstone. . . . . . . . . . . . . . . . . . . . . . .| 125.00000|
| |5|limestone. . . . . . . . . . . . . . . . . . . . . . .| 70.00000|
| |6|shale. . . . . . . . . . . . . . . . . . . . . . . . .| 125.00000|
| |7|sandy shale. . . . . . . . . . . . . . . . . . . . . .| 29.00000|
| |8|claysand . . . . . . . . . . . . . . . . . . . . . . .| 14.00000|
| |9|sand . . . . . . . . . . . . . . . . . . . . . . . . .| 272.00000|
|----------------------------------------------------------|----------|
|2|commercial and services | 115.00000|
| |--------------------------------------------------------|----------|
| |1|metamorphic. . . . . . . . . . . . . . . . . . . . . .| 16.00000|
| |4|sandstone. . . . . . . . . . . . . . . . . . . . . . .| 19.00000|
[...]
|----------------------------------------------------------|----------|
|8|transportation and utilities | 400.00000|
| |--------------------------------------------------------|----------|
| |4|sandstone. . . . . . . . . . . . . . . . . . . . . . .| 34.00000|
| |5|limestone. . . . . . . . . . . . . . . . . . . . . . .| 8.00000|
| |6|shale. . . . . . . . . . . . . . . . . . . . . . . . .| 104.00000|
| |7|sandy shale. . . . . . . . . . . . . . . . . . . . . .| 26.00000|
| |8|claysand . . . . . . . . . . . . . . . . . . . . . . .| 4.00000|
| |9|sand . . . . . . . . . . . . . . . . . . . . . . . . .| 224.00000|
|---------------------------------------------------------------------|
|TOTAL |1519.00000|
+---------------------------------------------------------------------+
模块可以显示一元统计值。它基于当前范围和栅格分辨率，计算像素数目和地图属性的最小值、最大值、平均值、中值、方差、标准偏差和方差系数：
g.region rast=elevation.10m -p
r.univar -g elevation.10m
n=2654802
min=1061.06
max=1846.74
range=785.679
mean=1348.37
stddev=175.494
variance=30798.3
coeff_var=13.0153
11.6. 处理栅格地图的方法¶
[klassen]
11.6.1. 重分类¶
[reclass]
在对栅格地图分类的时候，会创建一个新的属性表，原来的地图不受影响。内存的需求很小，因为这只涉及一个表。但原始地图对于新的分类地图来说是必需的基础图。
模块可以交互式操作，也可以通过命令行调用。必需的分类标准需要存储在一个文件中，该文件在模块执行时需要指明。
确切的命令语法应该在详细描述()中查找。例如，假定Spearfish中的道路需要重分类。我们不需要已有的五个分类，而只需指明状况的”好”或是”坏”：
# CLASSIFICATION BEFORE:
r.report roads
|1|interstate
|2|primary highway, hard surface
|3|secondary highway, hard surface
|4|light-duty road, improved surface
|5|unimproved road
|*|no data
# RECLASSIFICATION:
r.reclass in=roads out=roads.rcl
Enter rule(s), "end" when done, "help" if you need it
Data range is 1 to 5
> 1 2 3 = 1 good condition
> 4 5 = 2 bad condition
> end
# CLASSIFICATION AFTER:
r.report roads.rcl
|1|good condition
|2|bad condition
|*|no data
模块能够创建”独立的”地图：
g.region rast=roads.rcl -p
r.mapcalc "newmap=roads.rcl"
11.6.2. 掩模¶
[rastmask]
在栅格操作前进行掩模操作非常有用。掩模的设置影响所有后续的栅格分析，就像当前设置的范围和分辨率一样。
基本上，只有名为(大写)的地图才会被栅格模块认为掩模。在里，像素值被赋为(空值)的地方不会进行任何分析操作，在计算中使用除此之外的区域。
可以用多种不同的方法来创建掩模。如果有合适的地图，那么只需要拷贝或重命名。
g.copy rast=Mask,MASK
g.rename rast=Mask,MASK
模块提供了另一种可能性。它能够提取像素值，作为掩模。(请参阅[sec:mapcalc])
您有土地利用图，并且您只想分析高程在1200米以上的区域。为了确定这一限制条件，现选取高程地图，您能够用将它转换为掩模。掩模的创建如下：
g.region rast=elevation.10m -p
r.mapcalc "mask1200=if(elevation.10m > 1200.0,1,null())"
g.copy rast=mask1200,MASK
只有地图中表示为1的区域才会被分析，只要在当前地图集中。显示地图来验证它：
d.rast elevation.10m
d.rast -o roads
在任何一个项目中，掩模在命令行中的有效性可以通过来确认。
要删除掩模，使用命令。然后，整个当前范围都会被分析：
如果创建的掩模想在以后再次使用，为使其无效，只要重命名就可以了。
11.7. 栅格数据的数字化¶
[rdigit]
GRASS通过模块，可以数字化点、线和面。因此，每个要素能分配一个类别值和标签：
r.digit
Please choose one of the following
A define an area
C define a circle
L define a line
Q quit (and create map)
以这种方式创建掩模十分简单，但达不到所需的精度。