基於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);
select * from
--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）