[Học WebGIS nâng cao] – Bài toán tìm đường với PostGIS+ pgRouting

0
7769

Tiếp nối serie trước về Webgis sử dụng GeoServer, PostGIS, OpenLayer, series này sẽ giới thiệu các bạn một bài toán khá phổ biến với webgis là bài toán tìm đường đi ngắn nhất.

Tiếp nối các công cụ của bài trước (GeoServer, PostGreSQL, PostGIS, OpenLayer) thì bài này chúng ta sẽ cần thêm một công cụ nữa, đó là pgRouting, một extension của PostGIS cung cấp cho chúng ta giải pháp tìm đường với các thuật toán khác nhau. Các bạn có thể tìm hiểu thêm về pgRouting ở đây: http://pgrouting.org/

Cài đặt

Với các bạn sử dụng PostGIS bản mới nhất hiện nay (PostGIS 2.2.2  và PostGreSQL 9.5) thì pgRouting đã tích hợp sẵn vào bộ cài đặt, các bạn sẽ không phải cài thêm nữa, còn những bản cũ hơn ví dụ như bảnPostGreSQL 9.1 và PostGIS 2.0 không tích hợp sẵn pgRouting thì chúng ta sẽ tìm kiếm bộ cài tại đây: http://postgis.net/windows_downloads/ ví dụ mình tìm thấy link cho bản PostGreSQL 9.1 tại đây: http://winnie.postgis.net/download/windows/pg91/buildbot/. Sau đó chọn vào bản pgrouting-pg91-binaries-2.0.0devw64.zip để tải về.

Sau khi tải về chúng ta sẽ có 1 file zip, giải nén ra chúng ta sẽ được các thư mục bin, lib, share và các file khác, copy toàn bộ các thư mục và các file đó, paste vào thư mục cài PostGreSQL, ví dụ của mình là: C:\Program Files\PostgreSQL\9.1.

Tiếp theo chúng ta tạo Extension cho pgRouting:

Vào pgAdmin, chọn SQL rồi gõ lệnh sau

CREATE EXTENSION pgrouting;


Nếu các bạn cài PostGIS mà đã tích hợp sẵn pgRouting thì không cần làm bước này.

Kiểm tra xem bạn đã cài đặt thành công hay chưa và phiên bản pgRouting là bao nhiêu, bạn gõ lệnh sau

SELECT pgr_version();

Vậy là xong phần cài đặt.

Chuẩn bị dữ liệu

Để có dữ liệu cho bài toán tìm đường, chúng ta phải có một layer làm route dạng polyline.

Các bạn có thể tự tạo dữ liệu cho mình hoặc có thể tải của OpenStreetMap tại đây: http://download.geofabrik.de/asia.html . Sau khi tải xong chúng ta sẽ có một thư mục gồm nhiều file shp. Chúng ta sẽ sử dụng file roads.shp làm dữ liệu tìm đường.

Dữ liệu của OpenStreetMap là bản đồ đường đi của cả Việt Nam ( và nhiều nơi khác) nên khi sử dụng sẽ load hơi nặng. Nếu các bạn chỉ cần làm dữ liệu demo thì nên cắt ra một vùng nhỏ hơn để làm dữ liệu như hình dưới. Hệ quy chiếu của data là EPSG: 4236 WGS 84.

Clip thì chúng ta có thể dùng tool clip của QGIS. Đầu tiên chúng ta phải chuẩn bị một layer dạng polygon vẽ vùng chọn để clip (các bạn có thể vào Layer > Create Layer rồi chọn một loại layer chúng ta muốn tạo) . Sau đó các bạn vào Processing > ToolBox để bật toolbox. Trong toolbox chúng ta gõ clip vào ô tìm kiếm để mở chức năng clip

Chọn layer để clip, layer dùng làm polygon sẽ clip theo và file shp lưu ra. Ấn Run để hoàn tất việc clip

Lưu ý: Data của OpenStreetMap sẽ không thể sử dụng được luôn để đưa vào tạo route vì pgRouting đòi hỏi mỗi đoạn giao cắt sẽ phải là một line riêng. Vì vậy chúng ta sẽ phải tiến hành phá khối các polyline thành các line. Để phá khối polyline chúng ta sử dụng công cụ Explore của QGIS như sau:

  1. Vào Processing > ToolBox để bật toolbox. Trong toolbox chúng ta gõ Explore  vào ô tìm kiếm để mở chức năng Explore .
  2. Chọn layer để phá khối, chọn đường dẫn file sau khi đã phá khối và ấn Run

3. File shp mới sẽ được tạo ra với dữ liệu đã được phá khối.

Tiếp theo chúng ta sẽ đưa layer này vào database postGIS của chúng ta. Cách đưa các bạn có thể sử dụng QGIS như hướng dẫn sau: https://ungdungmoi.edu.vn/nang-cao-hieu-qua-voi-qgis.html. Chúng ta sẽ đặt tên bảng dữ liệu là roads.

Quay trở lại pgAdmin, để cho layer roads của chúng ta có thể tìm đường được chúng ta phải tạo  Topology cho nó theo các bước sau:

  1. Mở pgAdmin và chọn vào database chứa bảng roads và chọn SQL .

 

2. Chúng ta thêm 2 trường vào bảng roads tạo ở bước trước như sau:

alter table public.roads add column source integer;
alter table public.roads add column target integer;

 

3. Tạo topology cho roads như sau

select pgr_createTopology('public.roads', 0.0001, 'geom', 'id');

Các tham số như sau: tên bảng , độ phân giải, tên trường lưu geometry, tên trường id.
Ở đây mình chọn độ phân giải là 0.0001, các bạn có thể để độ phân giải khác.

4. Ấn Run để chạy lệnh trên. Nếu thực hiện thành công sẽ có thông báo như hình dưới:

5. Sau khi chạy xong, pgRouting sẽ tạo cho chúng ta một bảng nữa có tên là roads_vertices_pgr

Vậy là đã xong phần tạo dữ liệu, để test chức năng chúng ta sẽ sử dụng QGIS như sau:

  1. Mở QGIS, vào Database > DB Manager > DB Manager
  2. Chắc chắn rằng bạn đã connect đến PostGIS, chọn vào database bạn đã connect đến, vào menu Database > SQL window gõ lệnh sau rồi ấn Execute (F5)
SELECT seq, id1 AS node, id2 AS edge, cost, geom
FROM pgr_dijkstra(
'SELECT id, source, target, st_length(geom) as cost FROM public.roads',
1, 3000, false, false
) as di
JOIN public.roads pt
ON di.id2 = pt.id ;

3. Chọn Load as new layer để xem kết quả

hiển thị như hình dưới

 

Vậy là chúng ta đã tạo route thành công.

Tạo layer Route trong Geoserver

Ở trên chúng ta đã biết cách tìm đường với hàm pgr_dijkstra với cú pháp như sau: pgr_dijkstra(sql, source, target, directed). Hàm này sẽ tìm kiếm dựa vào 2 trường target và source trong bảng roads, đây là id của các điểm nút tạo được trong roads_vertices_prg. Để việc tìm kiếm trực quan hơn, ta cần tạo một hàm có thể tìm kiếm dựa vào tọa độ điểm đầu và tọa độ điểm cuối, đặt tên là pgr_fromAtoB như dưới đây:

CREATE OR REPLACE FUNCTION pgr_fromAtoB(
IN tbl varchar,
IN x1 double precision,
IN y1 double precision,
IN x2 double precision,
IN y2 double precision,
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
-- Find nearest node
EXECUTE 'SELECT id::integer FROM roads_vertices_pgr
ORDER BY the_geom <-> ST_GeometryFromText(''POINT('
|| x1 || ' ' || y1 || ')'',4326) LIMIT 1' INTO rec;
source := rec.id;
EXECUTE 'SELECT id::integer FROM roads_vertices_pgr
ORDER BY the_geom <-> ST_GeometryFromText(''POINT('
|| x2 || ' ' || y2 || ')'',4326) LIMIT 1' INTO rec;
target := rec.id;
-- Shortest path query (TODO: limit extent by BBOX)
seq := 0;
sql := 'SELECT id, geom, name, cost, source, target,
ST_Reverse(geom) AS flip_geom FROM ' ||
'pgr_dijkstra(''SELECT id , source::int, target::int, '
|| 'st_length(geom) as cost FROM '
|| quote_ident(tbl) || ''', '
|| source || ', ' || target
|| ' , false, false), '
|| quote_ident(tbl) || ' WHERE id2 = id 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.id;
name := rec.name;
cost := rec.cost;
geom := rec.geom;
RETURN NEXT;
END LOOP;
RETURN;
END;
$BODY$
LANGUAGE 'plpgsql' VOLATILE STRICT;

 

Hàm trên sẽ tìm trong roads_vertices_pgr để tìm các điểm gần nhất với điểm chúng ta nhập vào để làm target, source, rồi sử dụng 2 điểm đó đưa vào hàm pgr_dijkstra  để tìm route.

Hàm trên mình tham khảo trong hướng dẫn của pgRouting tại http://workshop.pgrouting.org/chapters/wrapper.html. Tuy nhiên có chỉnh sửa lại cho đúng với tên bảng, tên các trường dữ liệu, thay trường length trong tài liệu = st_length(geom) as cost. v..v. Các bạn cũng sẽ phải thay tên bảng roads_vertices_pgr, tên trường id, geom cho đúng với csdl của bạn.

Chúng ta copy hàm trên vào ô nhập sql của csdl, chỉnh sửa cho đúng với csdl của mình rồi ấn run, hàm này sẽ được nạp vào trong csdl.

Tiếp theo chúng ta vào trang quản trị GeoServer. Chúng ta tạo workspace và store cho layer theo hướng dẫn Public Data với GeoServer

Đến bước tạo layer, chúng ta sẽ dử dụng công cụ Configure new SQL view… của GeoServer. Chúng ta chọn Configure new SQL view…

Bạn điền View Name là route, SQL statement nhập như sau:

SELECT (route.geom) FROM (
SELECT geom FROM pgr_fromAtoB('roads', %x1%, %y1%, %x2%, %y2%
) ORDER BY seq) AS route

 

Tiếp theo chúng ta chọn Guess parameters from SQL, nhập Default value =0 và Validation regular expression = ^-?[\d.]+$

Trong phần Attributes, ấn Refesh, trong Type bạn chọn LineString, SRID chọn hệ tọa độ của data, ở đây là 4326.

Ấn Save để hoàn tất bước này. Tiếp theo chúng ta điền các thông tin khác cho layer. Bạn có thể không cần điền thêm gì, chỉ cần ấn Compute from data và Compute from native bounds để tạo tạo độ khung cho layer.

 

Ấn Save để hoàn tất việc tạo Layer.

Vậy là chúng ta đã tạo xong layer route trên GeoServer, bài sau chúng ta sẽ tìm hiểu cách sử dụng layer này để tìm đường với OpenLayer.

Tham khảo: http://workshop.pgrouting.org/chapters/geoserver.html

Lưu ý: hàm ST_MakeLine trong tài liệu của pgRouting không chạy được nên tạm thời mình bỏ đi, chúng ta có thể không cần gộp vào 1 line mà đưa ra nhiều geometry vẫn hiển thị được route.

Hiển thị route trên WebGIS với OpenLayer

Đầu tiên chúng ta sẽ sử dụng một layer roads để hiển thị đường làm bản đồ nền, các bạn đọc bài Ứng dụng WebGIS với Openlayer để xem cách tạo WebGIS sử dụng OpenLayer, GeoServer.

Layer của mình sẽ để như sau:

var format = 'image/png';
var bounds = [105.671539306641, 20.8914451599121,
105.982925415039, 21.186128616333];
roads = new ol.layer.Image({
source: new ol.source.ImageWMS({
url: 'http://localhost:8080/geoserver/dc/wms',
params: {
'FORMAT': format,
'VERSION': '1.1.1',
STYLES: '',
LAYERS: 'dc:roads',
}
})
});
var projection = new ol.proj.Projection({
code: 'EPSG:4326',
units: 'degrees',
axisOrientation: 'neu'
});
var view = new ol.View({
projection: projection
});
var map = new ol.Map({
target: 'map',
layers: [
roads
],
view: view
});
//map.getView().fitExtent(bounds, map.getSize());
map.getView().fit(bounds, map.getSize());

 

Chúng ta vẫn cần 1 thẻ div để hiển thị map

<div id="map" class="map"></div>

Chúng ta sẽ cần thêm 2 textbox để hiển thị tọa độ của điểm chúng ta chọn trên bản đồ và các nút tìm đường và xóa đường.

<input type="text" id="txtPoint1" />
<br />
<input type="text" id="txtPoint2" />
<br />
<button id="btnSolve">Tìm đường</button>
<button id="btnReset">Xóa đường</button>

Chúng ta sẽ viết thêm sự kiện click vào map để lấy tọa độ điểm chọn trên map và lưu vào biến startpoint và endpoint sau đó add 2 điểm này lên bản đồ bằng cách đưa vào 1 vector layer

var startPoint = new ol.Feature();
var destPoint = new ol.Feature();
var vectorLayer = new ol.layer.Vector({
source: new ol.source.Vector({
features: [startPoint, destPoint]
})
});
map.on('singleclick', function (evt) {
if (startPoint.getGeometry() == null) {
// First click.
startPoint.setGeometry(new ol.geom.Point(evt.coordinate));
$("#txtPoint1").val(evt.coordinate);
} else if (destPoint.getGeometry() == null) {
// Second click.
destPoint.setGeometry(new ol.geom.Point(evt.coordinate));
$("#txtPoint2").val(evt.coordinate);
}
});

Trong sự kiện click nút tìm đường, chúng ta sẽ gọi layer route chúng ta tạo ở trong GeoServer, đưa params là điểm đầu điểm cuối vào layer như sau:

var result;
$("#btnSolve").click(function () {
var startCoord = startPoint.getGeometry().getCoordinates();
var destCoord = destPoint.getGeometry().getCoordinates();
var params = {
LAYERS: 'dc:route',
FORMAT: 'image/png'
};
var viewparams = [
'x1:' + startCoord[0], 'y1:' + startCoord[1],
'x2:' + destCoord[0], 'y2:' + destCoord[1]
];
params.viewparams = viewparams.join(';');
result = new ol.layer.Image({
source: new ol.source.ImageWMS({
url: 'http://localhost:8080/geoserver/dc/wms',
params: params
})
});
map.addLayer(result);
});

Cuối cùng là sự kiện xóa đường như sau:

$("#btnReset").click(function () {
startPoint.setGeometry(null);
destPoint.setGeometry(null);
// Remove the result layer.
map.removeLayer(result);
});

Kết quả của chúng ta sẽ được như sau:

 

Code toàn bộ bài:

<html xmlns="http://www.w3.org/1999/xhtml">
<head >
<title>Openlayers test</title>
<link rel="stylesheet" href="http://openlayers.org/en/v3.15.1/css/ol.css" type="text/css" />
<script src="http://openlayers.org/en/v3.15.1/build/ol.js" type="text/javascript"></script>
<script src="https://code.jquery.com/jquery-1.12.3.min.js" type="text/javascript"></script>
<style>
.map, .righ-panel {
height: 500px;
width: 40%;
float: left;
}
.map {
border: 1px solid #000;
}
</style>
<script type="text/javascript">
var roads, result;
$("#document").ready(function () {
var startPoint = new ol.Feature();
var destPoint = new ol.Feature();
var format = 'image/png';
var bounds = [105.671539306641, 20.8914451599121,
105.982925415039, 21.186128616333];
roads = new ol.layer.Image({
source: new ol.source.ImageWMS({
url: 'http://localhost:8080/geoserver/dc/wms',
params: {
'FORMAT': format,
'VERSION': '1.1.1',
STYLES: '',
LAYERS: 'dc:roads',
}
})
});
var projection = new ol.proj.Projection({
code: 'EPSG:4326',
units: 'degrees',
axisOrientation: 'neu'
});
var view = new ol.View({
projection: projection
});
var map = new ol.Map({
target: 'map',
layers: [
roads
],
view: view
});
//map.getView().fitExtent(bounds, map.getSize());
map.getView().fit(bounds, map.getSize());
var vectorLayer = new ol.layer.Vector({
source: new ol.source.Vector({
features: [startPoint, destPoint]
})
});
map.addLayer(vectorLayer);
map.on('singleclick', function (evt) {
if (startPoint.getGeometry() == null) {
// First click.
startPoint.setGeometry(new ol.geom.Point(evt.coordinate));
$("#txtPoint1").val(evt.coordinate);
} else if (destPoint.getGeometry() == null) {
// Second click.
destPoint.setGeometry(new ol.geom.Point(evt.coordinate));
$("#txtPoint2").val(evt.coordinate);
}
});
$("#btnSolve").click(function () {
var startCoord = startPoint.getGeometry().getCoordinates();
var destCoord = destPoint.getGeometry().getCoordinates();
var params = {
LAYERS: 'dc:route',
FORMAT: 'image/png'
};
var viewparams = [
'x1:' + startCoord[0], 'y1:' + startCoord[1],
'x2:' + destCoord[0], 'y2:' + destCoord[1]
];
params.viewparams = viewparams.join(';');
result = new ol.layer.Image({
source: new ol.source.ImageWMS({
url: 'http://localhost:8080/geoserver/dc/wms',
params: params
})
});
map.addLayer(result);
});
$("#btnReset").click(function () {
startPoint.setGeometry(null);
destPoint.setGeometry(null);
// Remove the result layer.
map.removeLayer(result);
});
});
</script>
</head>
<body>
<div id="map" class="map"></div>
<div class="righ-panel">
<input type="text" id="txtPoint1" />
<br />
<input type="text" id="txtPoint2" />
<br />
<button id="btnSolve">Tìm đường</button>
<button id="btnReset">Xóa đường</button>
</div>
</body>
</html>

 

Chú ý: Bạn phải thiết lập style cho layer route để phân biệt với layer roads. Cách thiết lập style cho layer trong GeoServer các bạn xem ở đây nhé: Nâng cao hiệu quả với QGIS

Tham khảo: http://workshop.pgrouting.org/chapters/ol3_client.html

 

 

Tác giả: Đỗ Xuân Cường

Nguồn bài viết: cuongdx313.wordpress.com

Previous article[Xây dựng WebGIS #11] – Hiển thị Popup
Next articleỨng dụng hệ thống thông tin địa lý (GIS) – quản lý kinh tế xã hội

RELATED ARTICLESMORE FROM AUTHOR

LEAVE A REPLY

Please enter your comment!
Please enter your name here

Website này sử dụng Akismet để hạn chế spam. Tìm hiểu bình luận của bạn được duyệt như thế nào.