基于Arcgis Engine 10.2（C#）+PostgreSQL 11（Postgis 3）+pgRouting 3.0实现使用数据库进行路径规划
前言：最近在（被迫）使用ArcGIS Engine10.2（.NET平台）进行二次开发（桌面应用），因为想做一个最短路径查询的功能，而arcgis的网络分析又比较麻烦，于是想到了使用Postgis。但这样就需要将本地shp存到数据库，在程序中连接数据库。
百度了半天发现Arcgis Engine直接连接PostgreSQL数据库需要用到ArcSDE。ArcSDE还需要另外安装，而且我用的ArcGIS Engine10.2只支持PostgreSQL 9.x（我的数据库版本是11），这样似乎就很麻烦了啊（。）
行叭，另辟蹊径，因为Arcgis和Postgis肯定都遵循OGC的SFS规范，只需要在postgis中把最短路径规划函数写好，导出结果的wkb/wkt，然后用C#连接数据库（使用Npgsql）获取wkb/wkt，最后用ArcGIS Engine把wkb/wkt转成IGeometry类对象，显示到地图上就行了。
整体思路就是这样，其实很类似于WebGIS中路径规划的思路，下面详细介绍一下：
一、PostgreSQL的最短路径函数
这一步网上有很多文章，写的也还不错，但是大部分比较古老了，用的还是pgRouting 1.x或者2.0的函数，这给我造成了很多困扰。
这里我使用PostgreSQL11+postgis 3+pgRouting 3.0来实现。
1、准备数据
准备路网数据，需要将所有线段连接的地方全部断开，这样便于以后的分析。
如果数据在交叉处没有断开，可以在ArcgisDesktop中的ArcToolBox 找到数据管理工具（Data management tools）——>要素（Features）——>要素转线（Feature to line）
2、创建空间数据库，开启扩展
在新的空间数据库里添加空间扩展，使用如下SQL语句：
CREATE EXTENSION postgis;
CREATE EXTENSION pgrouting;
CREATE EXTENSION postgis_topology;
CREATE EXTENSION fuzzystrmatch;
CREATE EXTENSION postgis_tiger_geocoder;
CREATE EXTENSION address_standardizer;
3、导入数据
这一步应该都懂，就不过多介绍了。
值得注意的是，有人说导入时应该勾选下图的选项（生成简单要素），否则自动生成为Multi要素（多面、多线、多点）路径规划会失败。但是我自己做的时候没遇到这种问题，以防万一，可以勾上。
4、添加路网拓扑结构
这里我就只考虑无向图的情况，有向图是类似的可以自行百度。
--添加起点id
ALTER TABLE public.whu_road ADD COLUMN source integer;
--添加终点id
ALTER TABLE public.whu_road ADD COLUMN target integer;
--添加道路权重值
ALTER TABLE public.whu_road ADD COLUMN length double precision;
--为sampledata表创建拓扑布局，即为source和target字段赋值
SELECT pgr_createTopology('public.whu_road',0.0001, 'geom', 'gid');
--为source和target字段创建索引
CREATE INDEX source_idx ON whu_road("source");
CREATE INDEX target_idx ON whu_road("target");
--为length赋值
update whu_road set length =st_length(geom);
--为road_xblk表添加reverse_cost字段并用length的值赋值
ALTER TABLE whu_road ADD COLUMN reverse_cost double precision;
UPDATE whu_road SET reverse_cost =length;
5、创建最短路径函数，返回路径线段序列
其中最短路径使用dijkstra算法，调用pgrouting的pgr_dijkstra函数。
网上有很多教程使用pgr_kdijkstraPath，这种就是比较老的版本，pgrouting3.0是不能使用的。
--删除已存在的函数
DROP FUNCTION pgr_fromAtoB(tbl varchar,startx float, starty float,endx float,endy float);
--DROP FUNCTION pgr_fromAtoB(varchar, double precision, double precision
-- double precision, double precision);
--基于任意两点之间的最短路径分析
CREATE OR REPLACE FUNCTION pgr_fromAtoB(
IN tbl varchar,--数据库表名
IN x1 double precision,--起点x坐标
IN y1 double precision,--起点y坐标
IN x2 double precision,--终点x坐标
IN y2 double precision,--终点y坐标
OUT seq integer,--道路序号
OUT gid integer,
OUT name text,--道路名
OUT heading double precision,
OUT cost double precision,--消耗
OUT geom geometry--道路几何集合
)
RETURNS SETOF record AS
$BODY$
DECLARE
sql text;
rec record;
source integer;
target integer;
point integer;
BEGIN
-- 查询距离出发点最近的道路节点
EXECUTE 'SELECT id::integer FROM '|| quote_ident(tbl) ||'_vertices_pgr
ORDER BY the_geom <-> ST_GeometryFromText(''POINT('
|| x1 || ' ' || y1 || ')'',4326) LIMIT 1' INTO rec;
source := rec.id;
-- 查询距离目的地最近的道路节点
EXECUTE 'SELECT id::integer FROM '|| quote_ident(tbl) ||'_vertices_pgr
ORDER BY the_geom <-> ST_GeometryFromText(''POINT('
|| x2 || ' ' || y2 || ')'',4326) LIMIT 1' INTO rec;
target := rec.id;
-- 最短路径查询
seq := 0;
sql := 'SELECT gid, geom, node as name, cost, source, target,
ST_Reverse(geom) AS flip_geom FROM ' ||
'pgr_dijkstra(''SELECT gid as id,
source::integer,target::integer,'
|| 'length::float AS cost FROM '
|| quote_ident(tbl) || ''', '
|| source || ', ' || target
|| ' ,false) as di, '
|| quote_ident(tbl) || ' WHERE di.edge = gid ORDER BY seq';
-- Remember start point
point := source;
FOR rec IN EXECUTE sql
LOOP
-- Flip geometry (if required)
IF ( point != rec.source ) THEN
rec.geom := rec.flip_geom;
point := rec.source;
ELSE
point := rec.target;
END IF;
-- Calculate heading (simplified)
EXECUTE 'SELECT degrees( ST_Azimuth(
ST_StartPoint(''' || rec.geom::text || '''),
ST_EndPoint(''' || rec.geom::text || ''') ) )'
INTO heading;
-- Return record
seq := seq + 1;
gid := rec.gid;
name := rec.name;
cost := rec.cost;
geom := rec.geom;
RETURN NEXT;
END LOOP;
RETURN;
END;
$BODY$
LANGUAGE 'plpgsql' VOLATILE STRICT;
6、测试
select * from sqlpgr_fromAtoB('whu_road',114.3560,30.5362,114.3683,30.5484);
随便选两个点测试一下，OK，已经返回了路径的线段序列。
7、 路线合并，转换为wkb或wkt
因为最终我们需要将最短路径传到桌面端显示，为了方便，使用ST_Union将这些线段合并成一个完整的线路。使用方法如下：
select ST_Union(geom) as route from pgr_fromAtoB('whu_road'::text,114.353,30.539,114.367,30.544) ;
合并的路径我们看不出来是什么，不过将他转成wkt就能理解了：
select ST_astext(ST_Union(geom)) as route from pgr_fromAtoB('whu_road'::text,114.353,30.539,114.367,30.544) ;
当然，为了方便C#读取，还是使用ST_asbinary转换为wkb（二进制）比较好，这也是我们后面需要用到的语句：
select ST_asbinary(ST_Union(geom)) as route from pgr_fromAtoB('whu_road'::text,114.353,30.539,114.367,30.544) ;
二、C#连接读取数据库
到这里事情就简单了。
1、安装Npgsql
在VS中打开项目-管理NuGet程序包
搜索安装Npgsql
2、读取数据库
创建一个DAO类，写一个路经查询函数，这里给一个简单的例子：
namespace WHUGIS.Classes
{
class DAO
{
private static string connectionString = "User ID=postgres;Password=admin;Server=localhost;Port=5432;Database=GIS_engine;";
public DAO()
{
}
public static Byte[] executeRouteQuery(string sqlstr)
{
NpgsqlConnection sqlConn = new NpgsqlConnection(connectionString);
try
{
sqlConn.Open();
NpgsqlCommand objCommand = new NpgsqlCommand(sqlstr, sqlConn);
Byte[] routeWKB = (byte[])objCommand.ExecuteScalar();
return routeWKB;
}
catch(Exception ee)
{
MessageBox.Show(ee.Message);
return null;
}
finally
{
sqlConn.Close();
}
}
}
}
三、使用ArcGIS Engine在AxMapControl中绘制最短路径
让我们随便创建一个RouteForm，设计一个界面用于获取点坐标。
关于C#的界面设计和交互的实现我在这里就不多说了，下面是在AxMapControl绘制最短路径的简单实现（在RouteForm类当中）：
private IElement pElement = null;
//路径起止点
private double startX = 0, startY = 0;
private double endX = 0, endY = 0;
//主界面地图控件的引用
private AxMapControl mMapControl;
//确定按钮被点击函数
private void button1_Click(object sender, EventArgs e)
{
//连接数据库，路径规划
string sqlstr = null;
if ((startX * startY * endX * endY) != 0)
sqlstr = "select ST_asbinary(ST_Union(geom)) as route from pgr_fromAtoB('whu_road_feature2line'::text," + startX + "," + startY + "," + endX + "," + endY + ") ;";
else
{
MessageBox.Show("起始点或终止点不能为空");
return;
}
Byte[] routeWKB = DAO.executeRouteQuery(sqlstr);
IGeometry geom;
int countin = routeWKB.GetLength(0);
//地图容器，创建临时元素
IMap pMap = mMapControl.Map; //
IActiveView pActiveView = pMap as IActiveView;
IGraphicsContainer pGraphicsContainer = pMap as IGraphicsContainer;
if (pElement != null)
{
pGraphicsContainer.DeleteElement(pElement);
}
//转换wkb为IGeometry
IGeometryFactory3 factory = new GeometryEnvironment() as IGeometryFactory3;
factory.CreateGeometryFromWkbVariant(routeWKB, out geom, out countin);
IPolyline pLine = (IPolyline)geom;
//定义要素symbol
ISimpleLineSymbol pLineSym = new SimpleLineSymbol();
IRgbColor pColor = new RgbColor();
pColor.Red = 11;
pColor.Green = 120;
pColor.Blue = 233;
pLineSym.Color = pColor;
pLineSym.Style = esriSimpleLineStyle.esriSLSSolid;
pLineSym.Width = 2;
//线元素symbol绑定
ILineElement pLineElement = new LineElementClass();
pLineElement.Symbol = pLineSym;
//添加geom
pElement = pLineElement as IElement;
pElement.Geometry = pLine;
//加入地图并刷新
pGraphicsContainer.AddElement(pElement, 0);
pActiveView.Refresh();
}
最后让我们看看效果：
参考资料：
关于在ArcGIS平台使用PostGIS，对Geometry字段的处理
GeoServer+PostgreSQL+PostGIS+pgRouting实现最短路径查询
GeoServer+PostgreSQL+PostGIS+pgRouting实现最短路径查询
postgresql+postgis+pgrouting实现最短路径查询（1）－－－线数据的处理和建立拓扑...
PostgreSQL+PostGIS+pgRouting最短路径规划之
[AE] ArcGIS Engine - 地图整饰 - 添加图形元素（临时的Geometry）
