地理资源分析支持系统 (GRASS):不仅仅是一个地图绘制工具

地理资源分析支持系统 (GRASS):不仅仅是一个地图绘制工具

2016-10-27 作者: 肖芳 浏览: 2506 次

摘要: 当你下载安装包之后里面含 Oracle 支持的 GRASS 从其源代码编译 GRASS 需要安装多个支持库。由于大多数支持库均为各种 Linux 发布版本的标准(要获取强制性要求和可选要求的完整列表,请访问 http://grass.itc.it/grass6...

当你下载安装包之后里面含 Oracle 支持的 GRASS

从其源代码编译 GRASS 需要安装多个支持库。由于大多数支持库均为各种 Linux 发布版本的标准(要获取强制性要求和可选要求的完整列表,请访问 http://grass.itc.it/grass60/source/REQUIREMENTS.html),因此通常只需要单独下载和安装地图测绘投影库 (PROJ) 以及地理空间数据抽象库 (GDAL)/OGR 简单特性库。可以对 PROJ 使用预编译的二进制文件,而 GDAL/OGR 则需要从源代码进行编译才能包含 Oracle 支持。GDAL/OGR 是一个转换程序库能用于处理地理空间数据使用的大多数光栅和矢量格式。GRASS 将 GDAL/OGR 用作导入或导出后端,以最大限度地提高互操作性。

要从 GRASS 访问和写入 Oracle 数据库中存储的数据集,需要使用 Oracle 调用接口 (OCI)(默认情况下,该接口在 OGR 中不编译)编译 OGR。编译 GDAL/OGR 时,您需要定义到完整 Oracle 客户端库(标准 Oracle 安装提供了该库)的路径。由于某些头文件也需要编译,因此请确保下载完整的客户端库。

下载和编译包含 Oracle 支持的 GDAL/OGR 非常简单:

wget ftp://ftp.gdal.org/gdal/gdal-1.3.1.tar.gz 
tar xvfz gdal-1.3.1.tar.gz cd gdal-1.3.1 
./configure

--with-oci=/path/to/oci-library 
make 
make install

如果编译成功,则可以使用 ogrinfo 命令列出支持的矢量格式,如下面的简化示例所示:

ogrinfo --formats 
Loaded OGR Format Drivers:
-> "ESRI Shapefile" (read/write) 
[...]
-> "SQLite" (read/write) 
-> "ODBC" (read/write) 
-> "OCI" (read/write) 
-> "OGDI" (readonly) 
-> "PostgreSQL" (read/write) 
[...]

现在,OCI 绑定已经编译为 OGR,您可以使用刚刚安装的 GDAL/OGR 库继续编译 GRASS。要紧跟最新的开发,建议您下载 GRASS 6.1(支持使用 Oracle 数据库中存储的数据集)的每周 CVS 源代码快照。请按如下所示编译它。(有关源代码中各种编译选项的详细信息,请参阅 INSTALL 文件,网址是 http://grass.itc.it/grass61/source/snapshot/):

wget http://grass.itc.it/grass61/source/snapshot/grass-6.1.cvs_src_snapshot_[date].tar.gz 
tar xvfz grass-6.1.cvs_src_snapshot_[date].tar.gz 
cd grass-6.1.cvs 
./configure \ 
[...] --with-gdal=/usr/local/bin/gdal-config \ 
--with-proj
make 
make install

现在已经编译了所有内容,但您仍需要通过定义相关的环境变量 ORACLE_HOME 来定义 Oracle 数据库的安装路径:

export ORACLE_HOME=/path/to/oracle-installation 如果计划每天使用 Oracle 数据库,则可以在环境中永久定义环境变量 ORACLE_HOME:

echo "export ORACLE_HOME=/path/to/oracle-installation" >> ~/.bashrc

现在,您可以使用以下代码列出所有支持空间的数据库实例层,以对 Oracle 绑定进行测试

ogrinfo -summary -ro "OCI:user/pw@db_instance"

此时,可以将 GRASS GIS 与 Oracle 数据库一起使用来进行地理空间数据处理和分析。

开始使用现有数据集

学习 GRASS 最简单的方法是从使用现有的现成数据集开始。GRASS 网站提供了多个这样的数据集。本文使用的是南达科塔州某个地区的数据集(称作 Spearfish)。首先,创建一个用于存储所有 GRASS 数据的子目录(可以将它命名为 grassdata),然后将 Spearfish 数据下载并解压缩到该目录中:

cd $HOME
mkdir grassdata
cd grassdata
wget http://grass.itc.it/sampledata/spearfish_grass60data.tar.gz
tar xvfz spearfish_grass60data.tar.gz

GRASS 数据库和命令的结构。解压缩 Spearfish 数据集后,可以看到 GRASS 数据库的结构。GRASS 数据存储在一个称作数据库(也称作 GISDBASE)的目录中;在本示例中,该目录是您刚才在上面创建的 grassdata 目录。在该数据库中,项目按位置(数据库的子目录)进行组织:已下载的数据集包含位置 spearfish60。一定要注意的是,每个位置由它的坐标系、地图投影和地理边界定义。

每个位置可能有多个“地图集”,用于将项目划分为不同的主题、子区域或工作区,以供单个小组成员处理。每个地图集包含用于保存栅格、矢量和属性数据的子目录,以及一个名为 WIND 的工作空间范围定义文件。用户在 GRASS 中工作时看不到这些子目录和文件。定义新位置时,GRASS 将自动创建一个名为 PERMANENT 的特殊地图集,用于存储默认空间范围、坐标系定义,并还可能用于存储小组的地理数据。

地图集这个概念为团队协作提供了一个方便的环境,即允许任何数量的用户使用他们的地图集同时在一个位置中工作。仅当用户被授予权限时才能从其他地图集读取地图,但用户始终无法写入其他地图集。由所有小组成员共享的地图应存储在 PERMANENT 地图集中。利用这个简单的模式可以对大型 GIS 项目进行轻松的管理,甚至还可以使用远程安装的网络文件系统(如 NFS)进行管理。

GRASS 包含 350 多个模块,因此熟悉一下它的命名惯例很有必要。前缀指示功能组(命令执行的操作类型),句点后的单词描述了命令的作用或它使用的数据类型。

启动 GRASS 并显示数据。现在,您准备启动 GRASS 6.1;根据本地安装的不同,您可以通过菜单启动它,也可以通过键入以下内容从终端窗口启动

grass61 将打开一个图形用户界面。在第一个域中输入数据库的路径(本示例中为 /home/user1/grassdata)。然后,选择 spearfish60 作为位置,并选择 user1 作为地图集。

然后,单击左下角的 Enter GRASS。在显示以下有关 GRASS 版本以及如何获得帮助的基本信息后:

 Welcome to GRASS 6.1.cvs (2006)
GRASS homepage: http://grass.itc.it/
This version running through:   Bash Shell (/bin/bash)
Help is available with the command: g.manual -i
See the license terms with: g.version -c
If required, restart the graphical user interface with:d.m &
When ready to quit enter:exit
GRASS 6.1.cvs (spearfish60):~ >

您将收到 GRASS shell 提示,并将打开 GUI 显示管理器 (d.m)

如果您知道到数据库、位置和地图集的路径,则可以直接从命令行启动 GRASS:

grass61 /home/user1/grassdata/spearfish60/user1 尽管 GRASS 可以与多种类型的 GUI一起使用,但为了简单起见,本文使用了命令行。

获得帮助

要获得所有可用命令的概述以及在 HTML 浏览器中访问它们的手册页面,请键入

g.manual -i 如果您知道要获得其使用帮助的命令的名称,则可以通过使用 help 参数显示命令用法的简短描述。例如:

d.rast -help 要了解在您的位置可以访问哪些栅格和矢量数据,使用了哪种坐标系、地图投影和单位,以及空间范围(工作区的坐标)和当前栅格分辨率类型如何,请键入

g.list rast 
g.list vect 
g.proj -w 
g.region -p

如果您不知道投影的概念以及所有坐标的含义,现在也不必为此而感到担心。对于本文使用的示例,您不必了解其中的含义;但如果您想要(或需要)对地理空间数据执行一些实际的操作,那么至少需要了解一些基础知识(请参阅 http://www.gdf-hannover.de/lit_html/grass60_v1.2_en/index.html 中的相关部分)。

显示数据

GRASS d.m 或 QGIS 为显示复杂的栅格和矢量数据集提供了一个使用方便的 GUI,但对于本文的介绍性示例而言,使用一个简单的命令行显示足以满足需要:

 set the current spatial extent to the extent of the elevation raster map 
g.region rast=elevation.dem -p


start a graphic-monitor and display raster digital elevation model (DEM), 
 and roads (vector line data) 
d.mon x0 
d.rast elevation.dem d.vect roads

您还可以使用 GRASS 可视化模块 NVIZ 来浏览在三维空间中投影的数据。要了解有关可视化 GUI 用法的详细信息,请在键入以下代码后,从打开的 NVIZ 窗口右上角的 help 菜单中选择 NVIZ help

nviz elevation=elevation.dem vector=roads

分析关系

地理空间分析是 GRASS GIS 的主要的强项之一,您可以通过计算基本的地形参数(如坡度和方位)以及从数字标高模型 (DEM) 中提取河流网络和盆地边界来开始研究其功能。这些参数可以在很多应用场合(尤其是土地管理)中使用,并为个人使用提供了很有帮助的信息,例如在计划徒步旅行或自行车旅行时,您可以了解要穿过的河流数量以及途中不同地段的坡度。本示例结合使用了地形分析和土地所有权数据来查找控制河流沿岸土地的土地所有者。此类信息通常用于水土保持、水质保护或水资源调查。

 compute and display slope and aspect (terrain orientation) maps

r.slope.aspect elevation=elevation.dem slope=myslope

aspect=myaspect

d.rast myslope

d.rast myaspect

extract and display stream network and basins

r.watershed elevation=elevation.dem basin=mybasins \
stream=mystreams threshold=1000

d.rast mystreams

d.rast mybasins

您可以通过以下方法为提取的河流创建一个更复杂的表示:将河流的栅格表示细化为一个单元格的宽度,将其转换为矢量直线并导出它们的某些属性(例如河段长度):

 convert the raster stream map to vector line representation 
 and create an attribute table to store stream length attributes 
r.thin in=mystreams out=mystreams_thinned 
r.to.vect -s in=mystreams_thinned out=streams_vector 
v.db.addtable map=streams_vector table=streams_vector col="length int" 
d.vect streams_vector width=2

 create a new vector map with new categories added (vector IDs) 
 in this step an ID and table row is added to each vector line: 
v.category in=streams_vector out=streams_vector2 op=add 
v.to.db map=streams_vector2 op=cat column=cat

 verify that all IDs appear: 
v.db.select streams_vector2

 populate stream segment length to DB and display the length values on the map 
v.to.db map=streams_vector2 option=length column=length units=me 
v.db.select streams_vector2 
d.vect streams_vector2 
d.vect streams_vector2 display=attr attrcol=length

 optionally you can view the new basin and stream maps draped over terrain 
 in 3D space

nviz elevation.dem color=mybasins vect=streams_vector2

 clean up 
g.copy vect=streams_vector2,streams_final 
g.mremove vect="streams_vector*"

 store resulting stream-vectors in Oracle 
v.out.ogr in=streams_final format=OCI dsn="OCI:user/pw@db_instance" \
olayer=streams_oracle

 verification 
ogrinfo -summary -ro "OCI:user/pw@db_instance" streams_oracle
在下一步中,您可以结合使用河段与土地所有者数据,以找到河段沿岸土地的所有者并将每个土地所有者的河流总长度累加在一起。
 calculate intersection between landowners and stream-segments resulting in
 a new data set  streams_fields
v.overlay ain=streams_final atype=line bin=fields out=streams_fields operator=and

 add another column for segment-length per landowner
v.db.addcol streams_fields col="length int"

 add length to the newly created segments per fields
v.to.db map=streams_fields option=length column=length units=me
 add another column for accumulated segment-length per landowner
v.db.addcol streams_fields col="acclength int"

 export your intersected data set into Oracle
v.out.ogr in=streams_fields format=OCI dsn="OCI:user/pw@db_instance" \
olayer=streams_fields_oracle

GRASS 中最近新增的 OGR 绑定仍需要某些修复操作才能正确执行导出,因此您需要使用以下 SQL 语句手动调整特性 ID 列:

echo "UPDATE streams_fields_oracle SET \"OGR_FID\" = \"cat\";" \
| sqlplus user/pw@db_instance

现在,您可以将每个土地所有者的河段长度累加在一起。由于已经将数据集导出到 Oracle 中,因此可直接使用 sqlplus 命令行工具在 Oracle 中执行该操作。

 grab the minimum/maximum of category numbers from fields for updating the 
 segments
v.univar fields col=cat

 update the accumulated length for each landowner
for i in `seq 1 63`; do \
echo "UPDATE streams_fields_oracle SET \"acclength\" = \
(SELECT sum(\"length\") FROM streams_fields_oracle WHERE \"b_cat\"= $i) \
WHERE \"cat\" = (SELECT \"cat\" FROM streams_fields_oracle WHERE \"b_cat\"= $i \
AND rownum = 1);" >> /tmp/acclength.sql ; \
done

 update SQL in oracle
cat /tmp/acclength.sql | sqlplus user/pw@db_instance

 set the accumulated length to all related stream segments per landowner
for i in `seq 1 63`; do \
echo "UPDATE streams_fields_oracle SET \"acclength\" = (SELECT \
\"acclength\" FROM streams_fields_oracle WHERE \"b_cat\"=$i AND rownum = 1) \
WHERE \"acclength\" = 0 AND \"b_cat\" = $i;" >> /tmp/acclength_all.sql ;\
done

 update SQL for all stream segments
cat /tmp/acclength_all.sql | sqlplus user/pw@db_instance

更新 Oracle 数据库中的属性列后,可以使用 v.external 命令(将外部数据集“链接”到 GRASS 中)直接在 GRASS 中使用该数据。现在,您可以使用下列命令显示结果: re-link your modified data from Oracle back into GRASS v.external dsn=OCI:user/pw@db_instance out=streams_fields_external \ layer=streams_fields_oracle

 display the streams with over 2000m of accumulated length for a landowner
 in red with less than 2000m length in green along with the property
 boundaries (field) and basins
d.rast mybasins -o
d.vect fields type=boundary,centroid att=label lsize=10 xref=left \
yref=bottom display=shape,attr type=point,line,boundary,centroid size=1 \
layer=1 color=0:0:0 lcolor=0:0:0 fcolor=170:170:170 llayer=1 icon=basic/x 
d.vect streams_fields_external width=3 where="\"acclength\" >= 2000 AND \
\"b_label\" != 'Black Hills Natl.Forest' " col=red
d.vect streams_fields_external width=3 where="\"acclength\" < 2000 OR \
\"b_label\" = 'Black Hills Natl.Forest' " col=green

 display table
v.db.select streams_fields_external col=\"acclength\",\"b_label\" \
where="\"acclength\" >= 2000" |sort -nu

 cleanup
rm -f /tmp/acclength*

 exit GRASS

exit

结果显示,在 Spearfish 地区,最大的累计河流长度位于 Black Hills National Forest 内,为大多数河流提供了充分的保护。但有几个土地所有者对水质具有重要影响。要尝试了解有关这些土地所有者可能带来的影响的详细信息,您可以结合使用土地所有权地图字段和土地覆盖地图 landcover.30m 以及提取的河流,从而确定那些应鼓励其将河流沿岸的农用土地纳入水土保持计划的土地所有者。

如果将您自己的空间数据存储到 Oracle 数据库中,则您现在应能够读取该数据并将其导入到 GRASS 中,同时将地理空间分析的结果输出回 Oracle 中。您的数据很有可能位于与 spearfish60 位置 (UTM) 不同的坐标系中,因此您将必须了解如何在 GRASS 中创建一个新位置。文档“ Nutshell 中的 GRASS 6”提供了一个简要说明。

结束语 ##

本文从用户角度简要介绍了 GRASS GIS 系统,并演示了一个简单的空间分析。GRASS 是一款开放源代码的软件,其代码完全对外公开,因此可在其中修改模块或添加新模块。本代码用 C 语言编写,而如何在 GRASS 中进行编程则需要一篇单独的文章。如果您对开发很感兴趣,可以在 GRASS 网站中找到 GRASS 编程人员手册(PDF、HTML);该手册每周更新一次,且代码说明内嵌到了源代码文件(基于“doxygen”)中。您可以在 GRASS GIS 开发邮件列表中寻求帮助。

Helena Mitasova 是北加州大学海洋、地球和大气科学系的研究员,还是 GRASS 开发小组和开放源代码地理空间基金会的活跃成员。Markus Neteler 1999 年在德国汉诺威大学获得了自然地理和景观生态学硕士学位。他是 ITC-irst 的研究员以及意大利 Trento 的 CEA,并参与编写了两部有关 GRASS 的图书。

Stephan Holl 2003 年在德国汉诺威大学获得了自然地理和景观生态学硕士学位。他是德国 GDF Hannover 的项目经理,负责多个使用开放源代码 GIS 软件的项目。

关注“开源集思”公众号
获取免费资源

随机推荐


Copyright © from 2014. 开源地理空间基金会中文分会 吉ICP备05002032号

Powered by TorCMS

OSGeo 中国中心 邮件列表

问题讨论 : 要订阅或者退订列表,请点击 订阅

发言 : 请写信给: osgeo-china@lists.osgeo.org