前言：最近在（被迫）使用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）