基于ArcSDE的影像数据管理.doc

基于ArcSDE的影像数据管理疑惑篇 为了和解决篇组成一个完整的整体，决定两篇文章都以Article的方式发布。 文章转载请注明出处，毕竟是我的一篇技术文章 ESRI的ArcSDE是个性能比较强劲的空间数据引擎，它在管理矢量数据是采用的策略和方法效率都还是很不错的。我前一段做了一个基于ArcSDE的影像数据包装器，其作为一项WebGIS项目后端的一部分，功能就是完成前端对影像数据的任意范围查询，将通过ArcSDE存于Oracle Spatial中的影像数据实时提取出来并生成图片。基础影像数据的数据量为13GB。这里对影像数据做一下说明。影像数据一般有两种航拍照片和卫星照片，以航片为例（也就是本文所处理的对象），在入库前一般是若干张大小相同的已经数字化（也就是扫描好的）的带坐标信息的（好多定语啊）栅格图像文件，它们是彼此相邻的，是由飞机在一个区域进行“地毯式”拍照得到的。这种图像文件的数据量大约每个在几十兆左右，像素大小一般为几千点几千点。 在做的过程中我发现ArcSDE在管理影像数据方面就没有矢量数据那么强大了。我们项目的矢量部分是基于MO做的，但是MO明确不支持影像数据。所以我只能求助于SDE API（我们这里没有AO，不知道AO对影像支持的怎么样）。SDE API分成C-API和Java-API两种，其中又只有C-API支持影像（ESRI在搞什么飞机？），所以这就是我唯一的选择。 而使用C-API又是一件十分令人头疼的事情，因为文档非常有限，只有随软件的一份Developer Help，而这份文档用粗制滥造来形容一点都不冤枉它，其中信息量小，例子混乱，甚至还有错别字，实在是给ESRI砸牌子。这份文档在ESRI官方有下载，有兴趣的朋友可以去翻翻。 这份文档说明了SDE管理影像数据的方法，简而言之就是“两种方式，分层分块”。“两种方式”是指：栅格地图（Raster Map）和栅格目录（Raster Catalog）。 栅格地图比较适合带有标准坐标系的航片或卫片，当以这种方式导入影像文件的时候，SDE（实际上是ArcToolBox干的）会将它们拼接（mosaic）成一张完整的大地图，但这种方式对各个图幅的要求就非常严格，要实现拼接必须满足几个条件：相邻图幅的相邻边不能有一丁点重叠错位；图幅必须带坐标系等等。但这种把图幅全拼起来的管理方式也有问题，就是如果将来要更新其中的某一小块区域的图像怎么办呢？比如纽约的影像图，世贸大厦被撞没了，其实只需要更新一下曼哈顿那一点地方的数据就可以了，可要是以Raster Map方式管理全纽约的影像地图的话那么就意味着要全部删掉原来的数据，再重新导入一遍全城的航片，几十个G，不是小数目啊，如果是1:10000的话估计纽约的影像数据量应该在300个G左右，怎么也得导几天才能导完吧。 那么另一种方式呢？Raster Catalog说白了就是个相册。可以把风驴羊不相及的东西的照片都存进去，SDE只管存储和为它们维护一个目录，这一点可以从观察Raster Catalog对应的Oracle Spatial表看出来，感觉就像一个被肢解了的图片的集中营。当然ArcSDE毕竟是GIS软件，如果用它来存一些花花草草或者MM图片也确实划不来，因为它实在太慢了。但是如果存进去的是我前面说的相邻的若干航片，那么效果会有所不同。就体现在如果用ESRI的桌面工具（比如ArcMap，ArcCatalog）来预览这个Raster Catalog的时候，你将看到一张经过拼接的美图！我一直想知道这些工具背后实现拼接的算法，但未果。使用这种方式管理航片就没有前面Raster Map方式最后说的那种问题，比如世贸大厦没了，那么我们就重新在那里拍一张，再更新原来Raster Catalog中对应的那一个图幅就可以了，其它的则没必要去动。 前面说的“两种方式”是在比较高的尺度上讨论，在我所做的这个项目中，客户的基础数据采用了第二种方式，即Raster Catalog方式进行存储和管理。现在我们深入一点，进入图幅内部，看看它的存储和管理方式，这就是“分层分块”。 “分层”是指：影像金字塔（pyramid）索引。这个东西我不想多说，搞这方面的应该都知道，其基本思想就是利用采样自底向上生成金字塔，根据需求直接取其中某一级作为操作对象，以提高整体效率。当然就像这个世界中的其它事物一样，效率的提高是有代价的，这就是建塔带来的额外空间开销，建的级越多，越方便查询，当然数据冗余也越大。SDE可以为导入的影像建立金字塔索引。据我观察一般为6级。 “分块”是指每个图幅是按图块（block/tile）存储的，并使用格网索引。这里就又隐藏了一个巨大的问题，但还是听我说分块先。在将图幅存储于数据库中时，SDE不是傻乎乎的一行一行的存进去，而是将图幅划分成若干个大小相同的图块，每块大小不能超过16K，一般就取128128。划分的顺序是从上到下，从左至右的，分块的好处在于可以减少磁盘I/O。但是哪有那么好的图幅，长宽都是128的整数倍啊？而事实上就几乎没有这么好的图幅，那么图幅的长宽一般都是什么数呢？我手里的这些应该都是35002800的，但是它们的实际大小一般为356x285x，出现这种误差是很正常的，也是允许的。那这图幅的长和宽都除不开128怎么办呢？SDE的策略就是补零。在图幅的右侧和下侧补上若干零元（RGB:000000），也就是黑点，这样处理过的图幅就可以进行划分了，如图所示。 切下来的每一块作为一条记录再存入数据库。这样的策略导致问题随之而来：Raster Catalog中的每个图幅都有“黑边”，那又如何实现跨图幅的无缝拼接？前面我说了，ESRI的桌面工具可以做到，这就是我非常想知道它使用的算法的原因。那么我呢？目前唯一的希望就是SDE C-API能够提供读取这些“黑边”信息的接口，为此我给ESRI中国的技术支持发了E-mail。但结果是令人极度失望的，可爱的工程师在回信中热情的告诉我：C-API没有提供读取补零信息的接口，而且没有关于补零信息的任何文档。 至此，我的工作就停下了，我只能开始想别的出路。不过最终我手工实现了图幅去零和无缝拼接，这是后话了，具体方法将另文所述。 本文的目的就在于简单介绍一下ArcSDE管理影像数据的基本方法及其缺陷，在本文应用背景下，它的这种缺陷给我带来了很大的麻烦，虽然最终解决了问题，但还是感觉SDE自身应该可以采用更高明的影像管理方案，提供更方便的接口。单看影像这块，似乎不像ESRI这种公司做出来的东西，也许他压根就没想好好做这块。不知道同样领域有没有类似ArcSDE的替代品，更加好用，更方便做二次开发，还请各位指点，我指的是影像。 基于ArcSDE的影像数据管理解决篇 疑惑篇中简单介绍了基于ArcSDE的影像数据管理的基本方法、策略及其缺陷。那么要想基于ArcSDE的Raster Catalog实现对影像数据的任意范围查询，并且在跨图幅的情况下做到无缝拼接该怎么实现呢？ 我是这样做的。 首先说说逻辑上的思路。 问题的输入和输出都是很明确的。 输入：BBox, w 说明一下，BBox是Bounding Box的缩写，是以地理坐标表示的查询范围，由left, right, top和bottom四个参数组成；w为视口的大小，即最终显示在Web上（我说过的哦，这个影像数据源包装器是一个WebGIS项目服务器端的一部分）给用户看的地图窗口的大小，由width和height两个参数组成。 输出：与查询范围精确匹配的一幅w大小的jpeg图像文件 再说明一下，jpeg图像是经过有损压缩的，数据量比较小，我程序中每次得到的图像文件数据量一般在50100K左右，在Internet上传输是可以接受的。这种图像不带坐标信息，也无法进行二次处理，就是纯粹的一张图片，很纯的那种。这不太符合OGC的WCS规范，因为我现在只能提供jpeg图像；但却有点像WMS，如果把接口改改的话就算勉强符合。 那么根据ArcSDE中影像数据的管理和存储方式，从逻辑上我按照以下步骤完成功能： 第一步：根据BBox与w计算显示比例尺s； 第二步：根据s计算提取数据所在的金字塔级别l； 第三步：通过一次散列，求出与BBox相交的图幅序列ImageList(I1, I2, I3)； 第四步：将BBox分解到ImageList中的每一个图幅元素上成为子查询序列subBBox(sb1, sb2, sb3)； 第五步：在l级金字塔上，通过第二次散列，在Ii上求出与sbi相交的图块序列TileListi(B1, B2, B3)； 第六步：提取第三步得到的全部图块的数据； 第七步：将第四步得到的数据还原成位图文件 嗯，总体上就是这个样子。下面逐步细说。 第一、二步挺简单不用解释，第三步开始就有问题了。 在疑惑篇中我说了，SDE认为Raster Catalog就是一“相册”，其中存储的影像是没有什么关系的，更不会为它们建立什么索引结构。这样要想实现第三步中的散列就要先为Raster Catalog中的图幅建立索引表，这个表我建在了影像数据所在的表空间里，名字就叫IMAGEINDEX，里面为所有图幅按照它们之间的拓扑关系建立了格网索引每个图幅一个行号、一个列号。但这时就又有问题了，基础数据一共460个图幅，全拼起来并不是一个大长方形，而是基本与目标区域边界吻合的锯齿形。那么若建立格网索引就必然会有一些空的地方，就是说会有一些格网号并没有实际的图幅与之对应。那么对于这些空的地方，我是将它们也写入索引表，但标识其目标为0呢还是直接跳过它们不管呢？我选择了前者。这样，一个如下面示意图所示的格网索引就建起来了，第一次散列的问题也就可以解决了。而第五步中的第二次散列由于SDE对每个图幅中的图块是做了格网索引的（疑惑篇中有介绍），所以易于实现。 下面是最关键的第三到第五步，在这三步中，我要完成图幅的去零，以支持跨幅查询时的无缝拼接。因为我可以通过SDE API得到每个图幅原来的大小，即未补零时的大小，还可以得到图块的大小，通过它们就可以算出图幅边缘补零部分的宽度和高度，当查询涉及到那些带有零元的图块时，我就能知道其中有效信息部分的大小，进而只读取它们，跳过零元。但这带来一个小问题，就是图块的大小出现了不均匀。为此，我采用了一种比较笨的办法，就是让ImageList中的每个图幅Ii都自己记录自己的TileListi中每个图块的大小。这样在第六步提取数据和第七步恢复图像的时候就可以按照每个图块的实际大小来读写，而那些补的零就都可以去死了。这就是我去零的基本思路。 提取每个图块数据时，外层循环对ImageList的遍历顺序是Z序的，而内层循环在每个图幅内部对TileList的访问顺序也是Z序的。与之对应的，在恢复图像的时候，在每个图块内部绘制像素的顺序、在每个图幅内绘制图块的顺序、绘制图幅的顺序也都是Z序的。如下图所示： 最后一步得到位图还没完，需要再稍做加工，因为它的地理范围比BBox要大（因为提取数据的最小单位是图块而不是像素，这个应该好理解吧），所以要将它中间BBox对应的那一部分切出来。这还没完，切出来的部分可能不是w大小（这是因为金字塔索引是离散的，按比例尺拿数据只能就近取金字塔中的某一级），但不会差很远，所以还要稍微拉伸才行。这样两步加工之后得到的才是最终产品。 以上就是我整个思路的要点概况，其中关键的就在去零实现无缝拼接，不知道说清楚了没有。其实办法也挺初级的，疑惑篇的评论中bluntsword一下就给出解决思路了。本人比较笨，唉 下面再说说具体实现吧。事先说明一下，我的OOA/D经验相当匮乏，具体解决方案给出之后大家可能会觉得十分龌龊而不堪如目，其中可能会有不计其数的违背OO原则的地方，但请大家相信这也绝非我本意。诸位如果实在看不下去想骂两句的话就不用忍了。如果有一天我能看到某介绍设计方案的文献上有我的设计作为反面教材，那我也无怨无悔。不过虽然它很丑，但它在我们整个WebGIS系统中还算是比较稳定的一块。 SDE C-API是C写的（废话！），其中使用了少量C+的特性，不支持我相对比较熟悉的.net，所以我选择了VC6来实现。 这个影像数据包装器只负责从影像数据源中抽取数据，形成图片，并传回给GIS服务器，所以其实它连界面都不需要，但我建了个对话框应用，上面摆一个Edit，用于输出运行时的一些状态信息。 为实现功能，我建了一个CMyRaster类，它的内部完成了我想要的全部事情（这是不是一个典型的“全能类”啊，好怕怕）。其结构如下：（由于不能贴C+代码，所以这里暂且用C#格式的，下同） classCMyRaster public: staticWCHAR*ToWChar(char*str);/在GDI+中，有关字符的参数类型全部都是WCHAR类型，该函数是将传统字符串进行转换 staticintDisconnectSDE(SE_CONNECTION*connection);/断开与SDE的连接 staticintConnectSDE(char*SDE_servername,char*SDE_service,char*SDE_instance,char*SDE_user,char*SDE_password,SE_CONNECTION*connection);/连接SDE /这里应该提供一个ConnectSDE的重载版本用于连接SQLServer服务器 intGetRaster(constSE_CONNECTIONconnection,const_ConnectionPtrADOConn,doubleleft,doubleright,doubletop,doublebottom,longuserScale,CStringfilename);/提取影像数据并生成图片的核心方法 CMyRaster(); virtualCMyRaster(); private: intm_QueryRasterRight;/与查询范围相交的图幅序列的格网索引范围 intm_QueryRasterLeft; intm_QueryRasterTop; intm_QueryRasterBottom; RASTER_METADATA*m_RastersMetadata;/图幅元数据数组，核心变量 charstrTableNameSE_QUALIFIED_TABLE_NAME;/RasterCatalog的表名 charstrColNameSE_MAX_COLUMN_LEN;/数据的列名，一般就是IMAGE intm_QueryPyramidLevel;/查询所在的金字塔级别 SE_INTERPOLATION_TYPEm_Interpolation;/插值的方法 long*m_ScaleByLevel;/金字塔该级别上的比例尺 longm_NumOfBands;/波段数 longm_PyramidHeight;/金字塔高度 IMAGE_EXTENTm_DefaultTileSize;/tile大小的预设值 SE_ENVELOPEm_WholeExtent;/最大的全图范围 intm_UserMapWidth;/用户地图区的宽度 intm_UserMapHeight;/用户地图区的高度 intGet_RasterMetadata_by_BBox(constSE_CONNECTIONconnection,const_ConnectionPtrADOConn,doubleleft,doubleright,doubletop,doublebottom,longuserScale);/获取SDE中raster图层的元数据 intGet_RasterData_from_SDE(constSE_CONNECTIONconnection,const_ConnectionPtrADOConn);/从数据源中提取影像数据 intWrite_to_BMP_file(doubleleft,doubleright,doubletop,doublebottom,CStringfilename);/将影像数据恢复成位图文件 LONGGetScale(LFLOATras_xLength,LONGrasWidth,BOOLisBigFont);/计算比例尺 intGet_RasterMetadata_by_GridIndex(intrasrownbr,intrascolnbr,int*rasIndex);/根据格网索引取图幅 intGetImageCLSID(constWCHAR*format,CLSID*pCLSID);/得到格式为format的图像文件的编码值，访问该格式图像的COM组件的GUID值保存在pCLSID中 intGet_PyramidLevel_by_scale(longuserScale,int*pyramidLevel);/根据比例尺获取金字塔索引的级别 intGet_Result_extent(SE_ENVELOPE*resultExtent);/获取结果位图的地理范围 ; 从CMyRaster类的结构可以看出我将围绕 RASTER_METADATA *m_RastersMetadata; 这个动态数组来做文章。RASTER_METADATA是一个定义在CMyRaster外面的结构体，细节如下： typedefstructraster_metadata /下面是每个raster对应的唯一属性 intm_RasterID;/Raster_ID intm_RasterRowNum;/图幅格网索引行标号 intm_RasterColNum;/图幅格网索引列标号 intm_QueryTileLeft;/查询范围对应的tile最左一列的列号 intm_QueryTileRight;/查询范围对应的tile最右一列的列号 intm_QueryTileTop;/查询范围对应的tile最上一行的行号 intm_QueryTileBottom;/查询范围对应的tile最下一行的行号 SE_ENVELOPEm_RasterExtent;/全图的地理坐标范围 IMAGE_EXTENTm_ImageSize;/全图的像素范围 IMAGE_EXTENTm_QueryImageSize;/与查询范围相交的像素范围 IMAGE_EXTENT*m_QueryTilesSize;/查询范围对应的每个tile的大小（不带零） /下面是每个raster中每级金字塔对应的属性，均为数组，其大小在运行时才能确定 IMAGE_EXTENT*m_ImageSizeByLevel;/金字塔该级别上像素范围 IMAGE_EXTENT*m_ZerosSizeByLevel;/金字塔该级别上补零部分的宽与高 SE_ENVELOPE*m_RasterExtentByLevel;/金字塔该级别上的地理坐标范围 int*m_TotalTileColByLevel;/金字塔该级别上以tile为单位的列数 int*m_TotalTileRowByLevel;/金字塔该级别上以tile为单位的行数 RASTER_METADATA; 公有方法GetRaster()是完成功能的“主”方法，其中就这么三句话： /通过查询范围获取与查询范围相交的全部图幅的元数据 intret=Get_RasterMetadata_by_BBox(connection,ADOConn,left,right,top,bottom,scale); if(ret!=SE_SUCCESS) return-1; /获取上面求出的图幅序列的像素数据 ret=Get_RasterData_from_SDE(connection,ADOConn); if(ret!=SE_SUCCESS) return-1; /将存储在临时文件中的像素数据还原成位图 ret=Write_to_BMP_file(left,right,top,bottom,filename); if(ret!=SE_SUCCESS) return-1; 逻辑设计中的第一、二步都有与其对应的函数实现。第三至第五步则化为GetRaster()中的 /通过查询范围获取与查询范围相交的全部图幅的元数据 intret=Get_RasterMetadata_by_BBox(connection,ADOConn,left,right,top,bottom,scale); if(ret!=SE_SUCCESS) return-1; Get_RasterMetadata_by_BBox函数用于填充m_RastersMetadata数组，以供后面使用。 第六步提取数据由 /获取上面求出的图幅序列的像素数据 ret=Get_RasterData_from_SDE(connection,ADOConn); if(ret!=SE_