7. 栅格数据管理，查询和应用

7.1. 加载和创建栅格

大多数使用场景中，用户通常使用 raster2pgsql 将栅格文件导入 KingbaseES 中。

7.1.1. 使用 raster2pgsql 来加载栅格

加载程序 raster2pgsql 是一个可执行的栅格加载工具，它把支持 GDAL 格式的栅格文件转换成 SQL 化的文件，以便将其加载到 KGIS 栅格表中。

raster2pgsql 支持的栅格类型将和所依赖的 GDAL 库支持的类型是一样的。获取给定版本的 raster2pgsql 所支持的栅格类型，请使用 - G 参数。如果两边库一样的话，输出结果应该和 ST_GDALDrivers 输出结果一样。

注意

当使用函数 ST_CreateOverview 从一个已经对齐的栅格并指定参数 factor 因子创建 overview 时候，创建出的 overview 可以是不对齐的。

使用样例:

raster2pgsql raster_options_go_here raster_file someschema.sometable > out.sql

• -? 该参数展示帮助信息。如果执行 raster2pgsql 程序时候，不传递任何参数，也会展示帮助信息。

• -G 打印支持的栅格格式。

• (c|a|d|p) 下面这些选项互斥:

  • -c 创建新表，并把栅格数据导入到表中，这是默认模式。

  • -a 追加栅格数据到一个已经存在的表。

  • -d 先 drop 掉表，然后创建新表，并把栅格数据导入到表中。

  • -p 准备模式，只创建表，不导入数据。栅格处理过程：把适当的约束条件注册（写入）栅格系统表中。

• 栅格处理参数：为在栅格目录中正确注册所使用的约束：

  • -C 用于栅格约束注册 —— 比如栅格的 srid, pixelsize 元数据信息等等，保证栅格在视图 raster_columns 能够合适地注册。

  • -x 让设置最大边界约束失效，只有在 - C 参数同时使用时候才能使用这个参数。

  • -r 为规则的块设置约束（空间唯一性约束和覆盖瓦片）。只有在同时使用 - C 参数时，才能使用这个参数。

• 栅格处理参数：用于操作输入栅格数据集的可选参数：

  • -s 指定输出栅格使用指定的 SRID。如果没有提供或者提供值为 0，程序会检查栅格的元数据（比如栅格文件）以便决定一个合适的 SRID。

  • -b 从栅格要抽取的波段位置（下标从 1 开始）。想要抽取不止 1 个波段的话，使用逗号，分隔波段数字。如果没有指定，所有栅格的波段都会被抽取出来。

  • -t TILE_SIZE 瓦片大小。把每一个表的行表示的栅格切割成瓦片。TILE_SIZE 被表示成 WIDTH * HEIGHT 的形式，或者设置成值为 “auto” 来让加载程序使用第一个栅格参数计算合适的瓦片大小，然后把这个瓦片大小应用于所有的栅格。

  • -R, --register 把栅格以文件数据形式注册（栅格数据保存在数据库之外的文件系统中）。只有栅格的元数据和文件路径保存在数据库中（没有保存像素的数据）。

  • -l OVERVIEW_FACTOR 创建栅格的概览信息。如果 factor 值不只一个，使用逗号，进行分隔。创建的概览数据存储在数据库中，并且不受参数 - R 影响。注意生成的 SQL 文件会包含存储数据的主表和概览表。

  • -N NODATA 指定没有 NODATA 值的波段的 NODATA 值。

• 用于操作数据库对象的可选参数：

  • -q 把 KingbaseES 的 id 用引号包起来。

  • -f COLUMN 指定目标的栅格列名，默认是‘rast’。

  • -F 创建一个以文件名为列名的列。

  • -I 在栅格列上创建一个 GiST 索引。

  • -M 对栅格表执行 Vacuum analyze。

  • -T tablespace 指定新表的表空间。注意如果索引（包括 primary key）依然会使用默认的表空间，除非使用 - X 参数标识。

  • -X tablespace 指定表的索引表空间。如果使用了参数 - I，那么这个参数可以应用于 primary key 和空间索引。

  • -Y 使用 copy 语句而不是 insert 语句。

-e 单独地执行每一条语句，不使用一个大的事务执行所有语句。

-E ENDIAN 控制栅格输出时候生成的二进制的字节顺序，如果值为 0，表示取 XDR，如果值为 1，表示取 NDR（默认取值）。现在只支持 NDR 了。

-V version 指定输出版本的格式，默认是 0，当前值支持 0。

一个使用 raster2pgsql 的用于创建一个 sql 输入文件，并生成一个 100x100 方格的样例如下：

注意

可以省略表的 schema 名，比如 demelevation 而不是 public.demelevation ，这样栅格表就会创建在数据库默认的 schema 下面。

raster2pgsql -s 4326 -I -C -M *.tif -F -t 100x100 public.demelevation > elev.sqlpsql -d gisdb -f elev.sql

可以使用 UNIX 管道把转换和加载一步完成：

raster2pgsql -s 4326 -I -C -M *.tif -F -t 100x100 public.demelevation | ksql -d gisdb

这样就把栅格数据 Massachusetts 州的平面瓦片数据加载一个名为 aerials 的 schema 中，然后创建一个完整视图，然后使用 copy 语句把数据直接导入到数据库中。 参数 - e 不要求所有的插入语句都在一个事务中（这个参数的好处是可以让想立即在表中看到部分数据，因为分步提交可以不等待数据完全插入完成）。 程序会把栅格数据分成 128x128 的像素瓦片，然后加上栅格约束。这里使用 copy 而不是 insert 方式。加上参数 - F 会包含一个命名为从瓦片提取数据的文件名的字段。

raster2pgsql -I -C -e -Y -F -s 26986 -t 128x128 -l 2,4 bostonaerials2008*.jpg aerials.boston | psql -U kingbase -d gisdb -h localhost -p 5432

获取系统所支持的栅格类型列表：

raster2pgsql -G

参数 - G 输出如下列表：

Available GDAL raster formats:
Virtual Raster
GeoTIFF
National Imagery Transmission Format
Raster Product Format TOC format
ECRG TOC format
Erdas Imagine Images (.img)
CEOS SAR Image
CEOS Image
JAXA PALSAR Product Reader (Level 1.1/1.5)
Ground-based SAR Applications Testbed File Format (.gff)
ELAS
Arc/Info Binary Grid
Arc/Info ASCII Grid
GRASS ASCII Grid
SDTS Raster
DTED Elevation Raster
Portable Network Graphics
JPEG JFIF
In Memory Raster
Japanese DEM (.mem)
Graphics Interchange Format (.gif)
Graphics Interchange Format (.gif)
Envisat Image Format
Maptech BSB Nautical Charts
X11 PixMap Format
MS Windows Device Independent Bitmap
SPOT DIMAP
AirSAR Polarimetric Image
RadarSat 2 XML Product
PCIDSK Database File
PCRaster Raster File
ILWIS Raster Map
SGI Image File Format 1.0
SRTMHGT File Format
Leveller heightfield
Terragen heightfield
USGS Astrogeology ISIS cube (Version 3)
USGS Astrogeology ISIS cube (Version 2)
NASA Planetary Data System
EarthWatch .TIL
ERMapper .ers Labelled
NOAA Polar Orbiter Level 1b Data Set
FIT Image
GRIdded Binary (.grb)
Raster Matrix Format
EUMETSAT Archive native (.nat)
Idrisi Raster A.1
Intergraph Raster
Golden Software ASCII Grid (.grd)
Golden Software Binary Grid (.grd)
Golden Software 7 Binary Grid (.grd)
COSAR Annotated Binary Matrix (TerraSAR-X)
TerraSAR-X Product
DRDC COASP SAR Processor Raster
R Object Data Store
Portable Pixmap Format (netpbm)
USGS DOQ (Old Style)
USGS DOQ (New Style)
ENVI .hdr Labelled
ESRI .hdr Labelled
Generic Binary (.hdr Labelled)
PCI .aux Labelled
Vexcel MFF Raster
Vexcel MFF2 (HKV) Raster
Fuji BAS Scanner Image
GSC Geogrid
EOSAT FAST Format
VTP .bt (Binary Terrain) 1.3 Format
Erdas .LAN/.GIS
Convair PolGASP
Image Data and Analysis
NLAPS Data Format
Erdas Imagine Raw
DIPEx
FARSITE v.4 Landscape File (.lcp)
NOAA Vertical Datum .GTX
NADCON .los/.las Datum Grid Shift
NTv2 Datum Grid Shift
ACE2
Snow Data Assimilation System
Swedish Grid RIK (.rik)
USGS Optional ASCII DEM (and CDED)
GeoSoft Grid Exchange Format
Northwood Numeric Grid Format .grd/.tab
Northwood Classified Grid Format .grc/.tab
ARC Digitized Raster Graphics
Standard Raster Product (ASRP/USRP)
Magellan topo (.blx)
SAGA GIS Binary Grid (.sdat)
Kml Super Overlay
ASCII Gridded XYZ
HF2/HFZ heightfield raster
OziExplorer Image File
USGS LULC Composite Theme Grid
Arc/Info Export E00 GRID
ZMap Plus Grid
NOAA NGS Geoid Height Grids

7.1.2. 使用 KGIS 栅格函数来创建栅格

在许多场景下，会想在数据库中创建栅格即栅格表。KGIS 提供了相当多的函数来干这个事。一般步骤如下：

1. 创建一个带有栅格记录的栅格列的表可以用下面的完成：

   CREATE TABLE myrasters(rid serial primary key, rast raster);

2. 有许多函数可以完成这一步骤，如果创建的栅格不依赖于其他栅格，那么可以使用函数：ST_MakeEmptyRaster，接着使用 ST_AddBand 也可以使用 geometry 对象来创建栅格。 要想实现这个目的，需要使用函数 ST_AsRaster。可能还需要和其他函数比如函数 ST_Union 或函数 ST_MapAlgebraFct 或者其他地图代数系列函数联合使用。 甚至还有一些根据一些已经存在的栅格表创建新的栅格表的可选函数。例如可以使用函数 ST_Transform 根据一个已有的栅格表在其他投影系中创建一个新的栅格表。

3. 一旦完成了初始化栅格表，需要在栅格列上创建一个空间索引，使用如下语法：

   CREATE INDEX myrasters_rast_st_convexhull_idx ON myrasters USING gist( ST_ConvexHull( rast) );

   注意这里函数 ST_ConvexHull 的使用，因为大多数栅格操作符都是基于栅格的凸包。

4. 使用函数 AddRasterConstraints 添加约束。

7.2. 栅格系统目录

和 KGIS 打包在一起的有两个栅格视图。这两个视图都利用嵌入在栅格约束中的信息进行展示。因此这两个视图总是与表中的栅格数据保持一致性兼容，因为表中强加了约束。

1. raster_columns 这个视图记录了包含数据库中所有的栅格表的列。

2. raster_overviews 这个视图记录了包含的数据库中所有的栅格表的列，并以细粒度表的方式概述。这种类型的表在使用栅格加载程序并加上 - l 参数时候将会生成。

7.2.1. 栅格列的 Catalog（目录）

raster_columns 是一个包含的数据库所有的栅格表的栅格类型的列的 catalog 目录。这个视图利用嵌入在栅格约束中的信息进行展示。 因此这两个视图总是与表中的栅格数据保持一致性兼容，即便从一个其他数据库的备份中恢复栅格表。 如果在创建表时候没有使用加载工具，或者使用加载工具时候没有指定 - C 标识，可以使用函数 AddRasterConstraints 来强加约束， 这样的话 raster_columns 目录就可以将有关栅格瓦片的信息进行注册。

• r_table_catalog 表示数据表所在的库。它总会读取当前数据库信息。

• r_table_schema 表示栅格表所属的 schema。

• r_table_name 栅格表。

• r_raster_column 表名为 r_table_name 的栅格类型列。KGIS 中没有设置让每个表有多个栅格列，因此有可能一个栅格表会有多个不同的栅格列。

• srid 栅格表的 SRID，可以在 3.3.1 节了解相关。

• scale_x 空间参考坐标和像素的缩放比例。只有在栅格列的所有瓦片的 scale_x 是相同的才可用，然后会加上这个约束。请参考函数 ST_ScaleX 获取更多细节。

• scale_y 空间参考坐标和像素的缩放比例。只有在栅格列的所有瓦片的 scale_y 是相同的才可用，然后会加上这个约束。请参考函数 ST_ScaleY 获取更多细节。

• blocksize_x 每个栅格瓦片的宽度 (像素横向的个数)，参考函数 ST_Width 获取更多细节。

• blocksize_y 每个栅格瓦片的高度 (向下的像素数)。参考 ST_Height 获取更多细节。

• same_alignment 如果所有的栅格瓦片有相同的 alignment 对齐方式，那么返回 true。参考函数 ST_SameAlignment 获得更多细节。

• regular_blocking 如果栅格列有空间唯一和覆盖瓦片的约束，那么参数值应该是 TRUE，否则是 FALSE。

• num_bands 栅格集合的每个瓦片的波段数。这个与函数 ST_NumBands 返回结果一样。

• pixel_types 定义每一个波段的 pixel 类型的数组。这个数组元素的个数和波段的总数是一样的。这个参数是函数 ST_BandPixelType 的参数之一。

• nodata_values 表示每一个波段的 nodata_value 的 double 型数组。这个数组元素的个数和波段的总数是一样的。大多数操作中定义每个波段 pixel 值可以忽略。这个参数和函数 ST_BandNoDataValue 返回的信息类似。

• extent 这个是的栅格集合的所有栅格的边界。如果准备加载更多数据，但却会改变栅格的边界，那么需要运行函数 DropRasterConstraints 删除约束再加载数据，加载完成后，再运行函数 AddRasterConstraints 添加约束。

• spatial_index 如果一个栅格列有空间索引，那么返回值为 true。

7.2.2. 栅格概视图

raster_overviews 记录了用于概视图的栅格表的列的信息以及额外的何时需要利用概视图的有益信息。概视图表记录在 raster_columns 和 raster_overviews 中，因为这是本身他们的功能，但同样可以用于把低分辨率作为高分辨率的简化的目的。

当在加载数据，使用 - l 开关时候，他们会沿着主栅格表生成的。

概视图表包含了和其他栅格表一样的约束以及额外的对于概视图的信息。

注意

raster_overviews 里面的信息不是重复 raster_columns 里面的信息。 如果想得到 raster_columns 里面的的概视图的信息，可以把 raster_overviews 和 raster_columns 做关联，就可以获得全部信息。

使用概视图的两个主要理由是：

1. 中心表的分辨率的描述通常用于快速的地图放大。

2. 通常这种情况下会计算比更高分辨率的父瓦片计算更快一些，因为分辨率更低的会有更少的记录数，每个像素会覆盖更多的区域。虽然计算没有更高分辨率的精确，但是在许多经验规则下通常这种分辨率是足够的了。

raster_overviews 包含了如下列的信息：

• o_table_catalog 概视图表所在的数据库。它总是读取当前数据库信息。

• o_table_schema 栅格表概视图所属的数据库 schema。

• o_table_name 栅格概视图表名。

• o_raster_column 概视图表的栅格列。

• r_table_catalog 概视图服务的栅格表所在的数据库。它总是读取当前数据库信息。

• r_table_schema 概视图提供服务的栅格表所在的数据库 schema。

• r_table_name 概视图服务的栅格表名。

• r_raster_column 概视图提供信息的栅格列。

• overview_factor 这个是概视图表的金字塔级别。这个值越大，概视图表的分辨率越低。如果给定了一个图片的目录，加载工具 raster2pgsql 会计算每个图片文件的概视图，然后单独加载。 这个参数默认值是 Level 1，并且针对的是原始的文件。Level2 会让每个瓦片代表 4 个原始文件。 因此举个例子，如果有一个包含 5000x5000 像素图片的文件夹，并且打算分成每个 125x125 像素的块，因此对于的基础表的每个图像文件，将会得到 (5000*5000)/(125*125) = 1600 条记录， 其中的 (level=2) o_2 表将会有 ceiling (1600/Power (2,2)) = 400 条记录， 其中的 (level=3) o_3 表将会有 ceiling (1600/Power (2,3)) = 200 条记录。 如果像素数不能用瓦片的大小除尽，会得到一些碎片瓦片（瓦片没有完全填满）。 注意 raster2pgsql 生成的每个概览视图瓦片都和该瓦片的父级瓦片有相同的像素个数，但是分辨率更低，因为这个瓦片的每个像素相当于代表原始的 Power (2,overview_factor) 这么多个像素。

7.3. 创建定制应用

KGIS 栅格功能为提供了许多 SQL 函数用于把栅格与许多已知格式的数据进行交换，事实上这让有了许多数据交换的选择。例如可以使用办公工具 OpenOffice / LibreOffice 来展示栅格数据。除此之外，可以使用一些广泛使用的语言来演示本节的样例。

7.3.1. 在 PHP 中使用函数函数 ST_AsPNG 和其他栅格函数输出栅格的样例

本节，将会展示怎样使用 PHP KingbaseES 驱动和 ST_AsGDALRaster 系列函数来输出一个栅格的波段 1,2,3 到一个 PHP 请求数据流中，这样就可以把它嵌入到一个 HTML 标签 img 的 src 属性中。

在这个样例中，需要使用驱动 npgsql .NET KingbaseES，可以从这里获取到最新的驱动版本 http://npgsql.projects.KingbaseES。

只需要下载最新的版本并且把它丢进 ASP.NET 的 bin 目录就可以了。下面的查询样例展示了怎样联合使用栅格的一系列函数来获取到一个特殊的 wgs 84 缓冲区相交的所有瓦片，然后使用函数 ST_Union 对相交的瓦片求并集，然后把返回的所有波段使用函数 ST_Transform 转换到用户自定义的投影系中，最后再使用函数 ST_AsPNG 把它输出到一个 PNG 图片文件中。

会调用下面的 URL

http://mywebserver/test_raster.php?srid=2249

来获取 Massachusetts state plane feet（Massachusetts 州标准单位）的栅格图像

7.3.2. 在 ASP.NET C# 中使用函数函数 ST_AsPNG 和其他栅格函数输出栅格的样例

只需要下载最新的版本并且把它丢进 ASP.NET 的 bin 目录就可以了。下面的查询样例展示了怎样联合使用栅格的一系列函数来获取到一个特殊的 wgs 84 缓冲区相交的所有瓦片，然后使用函数 ST_Union 对相交的瓦片求并集，然后把返回的所有波段使用函数 ST_Transform 转换到用户自定义的投影系中，最后再使用函数 ST_AsPNG 把它输出到一个 PNG 图片文件中。

下面的样例和 6.3.1 节类似，只是使用了 C# 进行实现：

会调用下面的 URL

http://mywebserver/TestRaster.ashx?srid=2249

来获取 Massachusetts state plane feet（Massachusetts 州标准单位）的栅格图像

-- web.config connection string section --



// Code for TestRaster.ashx
<%@ WebHandler Language="C#" Class="TestRaster" %>
using System;
using System.Data;
using System.Web;
using Npgsql;
public class TestRaster : IHttpHandler
{
public void ProcessRequest(HttpContext context)
{
context.Response.ContentType = "image/png";
context.Response.BinaryWrite(GetResults(context));
}
public bool IsReusable {
get { return false; }
}
public byte[] GetResults(HttpContext context)
{
byte[] result = null;
NpgsqlCommand command;
string sql = null;
int input_srid = 26986;
try {
using (NpgsqlConnection conn = new NpgsqlConnection(System.Configuration. ←
ConfigurationManager.ConnectionStrings["DSN"].ConnectionString)) {
conn.Open();
if (context.Request["srid"] != null)
{
input_srid = Convert.ToInt32(context.Request["srid"]);
}
sql = @"SELECT ST_AsPNG(
ST_Transform(
ST_AddBand(
ST_Union(rast,1), ARRAY[ST_Union(rast,2),ST_Union(rast,3)])
,:input_srid) ) As new_rast
FROM aerials.boston
WHERE
ST_Intersects(rast,
ST_Transform(ST_MakeEnvelope(-71.1217, 42.227, ←
-71.1210, 42.218,4326),26986) )";
command = new NpgsqlCommand(sql, conn);
command.Parameters.Add(new NpgsqlParameter("input_srid", input_srid));
result = (byte[]) command.ExecuteScalar();
conn.Close();
}
}
catch (Exception ex)
{
result = null;
context.Response.Write(ex.Message.Trim());
}
return result;
}
}

7.3.3. 将栅格查询通过 Java 控制台程序输出成图像文件

这是一个简单的 Java 控制台程序，先做返回一个图像的查询，然后输出成指定文件。

可以通过如下命令 “编译” 自己的 Java 代码：

javac SaveQueryImage.java
jar cfm SaveQueryImage.jar Manifest.txt \*.class

使用命令调用这个函数方法如下：

java -jar SaveQueryImage.jar "SELECT ST_AsPNG(ST_AsRaster(ST_Buffer(ST_Point(1,5),10, ’ quad_segs=2’),150, 150, ’8BUI’,100));" "test.png"

-- Manifest.txt --
Class-Path: KingbaseES.jdbc4.jar
Main-Class: SaveQueryImage

// Code for SaveQueryImage.javaimport java.sql.Connection; import
java.sql.SQLException; import java.sql.PreparedStatement; import
java.sql.ResultSet;
import java.io.*;
public class SaveQueryImage {
public static void main(String[] argv) {
System.out.println("Checking if Driver is registered with DriverManager.");

try {
//java.sql.DriverManager.registerDriver (new
org.KingbaseES.Driver());Class.forName("org.KingbaseES.Driver");
}
catch (ClassNotFoundException cnfe) {
System.out.println("Couldn’t find the driver!");
cnfe.printStackTrace();
System.exit(1);
}

Connection conn = null;

try {
conn =
DriverManager.getConnection("jdbc:KingbaseES://localhost:5432/mydb","myuser
", "mypwd");
conn.setAutoCommit(false);

PreparedStatement sGetImg = conn.prepareStatement(argv[0]);

ResultSet rs = sGetImg.executeQuery();

FileOutputStream fout;
try{
rs.next();
/*\* Output to file name requested by user \**/fout = new
FileOutputStream(new File(argv[1]) );fout.write(rs.getBytes(1));
fout.close();
}
catch(Exception e)
{
System.out.println("Can’t create file");
e.printStackTrace();
}
rs.close(); sGetImg.close();
conn.close(); }
catch (SQLException se) {
System.out.println("Couldn’t connect: print out a stack trace and
exit.");
se.printStackTrace();
System.exit(1);
}
}
}

7.3.4. 通过 SQL 使用 PLPython 来导出 image 图像

这是一个 plpython 语言函数，它对每一行记录在服务器目录上创建一个文件。前提条件是：

• KingbaseES 服务器安装在 Linux 服务器上。

• KingbaseES 服务器安装了 plpython 扩展。这个函数可以在 plpythonu 和 plpython3u (Python 版本是 3.x) 下使用。

• 文件系统中存在一个名为 “/home/user/temp” 的目录，并且用户对该目录有读写权限。

CREATE OR REPLACE FUNCTION write_file (param_bytes bytea, param_filepath text)
RETURNS textAS $$
f = open(param_filepath, 'wb+')
f.write(param_bytes)
return param_filepath
$$ LANGUAGE plpythonu;

-- write out 5 images to the KingbaseES server in varying sizes
-- note the KingbaseES daemon account needs to have write access to folder
-- this echos back the file names created;
SELECT write_file(
   ST_AsPNG(ST_AsRaster(ST_Buffer(ST_Point(1,5),j*5, 'quad_segs=2'),150*j, 150*j, '8BUI',100)),
   '/home/user/temp/slices'|\| j \|\| '.png')
FROM generate_series(1,5) As j;
write_file
---------------------
/home/user/temp/slices1.png
/home/user/temp/slices2.png
/home/user/temp/slices3.png
/home/user/temp/slices4.png
/home/user/temp/slices5.png

7.3.5. 使用 PSQL 输出栅格

不幸的是 PSQL 没有简单的内置功能来输出二进制栅格数据。虽然有点 hack，但基于如下文章的建议 Clever Trick Challenge --Outputting bytea with psql，可以提供对大型对象导出的支持。要使用这个方法，首先通过命令行方式连接数据库。不像 Python 方法，这个方法在本地的电脑上创建文件。

SELECT oid, lowrite(lo_open(oid, 131072), png) As num_bytes
FROM
( VALUES (lo_create(0),
ST_AsPNG( (SELECT rast FROM aerials.boston WHERE rid=1) )
) ) As v(oid,png);

-- you’ll get an output something like --

oid       | num_bytes
----------+-----------
2630819   | 74860

-- next note the oid and do this replacing the c:/test.png to file path location
-- on your local computer
 \lo_export 2630819 ’C:/temp/aerial_samp.png’
-- this deletes the file from large object storage on dbSELECT lo_unlink(2630819);