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

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

Tham khảo bài gốc ở :
[Học WebGIS nâng cao] – Bài toán tìm đường với PostGIS+ pgRouting

Nhận xét

Đăng nhận xét

Bài đăng phổ biến từ blog này

DownLoad ảnh vệ tinh Landsat 8

Hình ảnh
Để bắt đầu với phần mềm Envi ta cần phải có dữ liệu vệ tinh để tiến hành xử lí. Như vậy, ta có thể Download ở đâu ???? Hiện nay, có rất nhiều website chia sẻ dữ liệu vệ tinh miễn phí như: +Website của cục khảo thí địa chất Hoa Kỳ (USGS): https://earthexplorer.usgs.gov/ +Website chia sẻ dữ liệu MODIS của UWG: https://ladsweb.modaps.eosdis.nasa.gov/ +Website chia sẻ dữ liệu Sentinel của ESA: https://scihub.copernicus.eu/ …………… Nếu muốn tìm hiểu sâu các bạn có thể truy cập vào một trong những web bên trên để có thể kiếm thêm thông tin về dữ liệu mà các bạn muốn sử dụng phục vụ nghiên cứu. Để làm ví dụ cơ bản ta sẽ sử dụng bộ dữ liệu ảnh quang Landsat 8 để tiến hành download (USGS) Đặc điểm ảnh Landsat 8: Là bộ dữ liệu mới nhất của Vệ tinh Landsat, với độ phân giải trung bình (30x30m) đã được tiền xử lí sơ bộ, mỗi tệp ảnh Landsat 8 được down về được đăng kí tọa độ WGS84 và hiệu chỉnh yếu tố địa hình rất thuận tiện để sử dụng. +Bắt đầu tải ảnh với (USGS), cần thực hiện đăng ký một
Đọc thêm

CẮt ảnh theo ranh giới (Sử dụng File ranh giới dạng *shp)

Hình ảnh
Công tác cắt ảnh này giúp loại bỏ những thông tin thừa ra khỏi khu vực nghiên cứu của bạn, Hiện nay trên youtube cũng có hướng dẫn về phần này rất nhiều. Tuy nhiên để cho các bạn hiểu hơn về quy trình cắt ảnh, mời cái bạn xem bài viết sau.   Đầu tiên: để cắt được ảnh ta cần vùng ranh giới dưới dạng *.shp.  Vì lí do các bài trước được viết dựa trên dữ liệu ảnh của khu vực TPHCM, cho nên mình sẽ lấy ranh giới TPHCM dưới dạng *.shp để thực hiện quá trình cắt ảnh. B1. Ranh giới dưới dạng *.shp khi mở bằng phần mềm ArcGis có dạng như sau B2. Mở Envi Classical 5.3, mở ảnh đã ghép sẵn chứa khu vực cần cắt. Trên thanh menu Chọn File/Open Vector File . Chọn đường dẫn vào Folder chứa thư mục *.shp của bạn. Rồi nhấn Open >>> Xem lại: Ghép ảnh vệ tinh B3. Hiển thị hộp thoại ta chọn Choose để chọn đường dẫn lưu file B4. Hiển thị ra hộp thoại chứa file *.evf ,  ta Click vô file vừa mới chuyển nhấn Load Selected/ Display 1#. Ta Click vô  File/Export Active Layer to ROIs
Đọc thêm

Bản đồ du lịch Việt Nam

Hình ảnh
Bản đồ du lịch Việt Nam thể hiện đầy đủ các tuyến điểm du lịch cần thiết với các thông tin bản đồ du lịch cho từng tỉnh/thành, khu vực riêng biệt. Đối với bản đồ du lịch tỉnh hiện được được cập nhật, hỗ trợ mới nên chưa thống nhất. Bản đồ du lịch tuyến điểm Việt Nam [caption id="attachment_1657" align="alignnone" width="2000"] Bản đồ du lịch Việt Nam dành cho dân phượt bụi[/caption] Bản đồ du lịch tuyến du lịch nổi tiếng Việt Nam [caption id="attachment_1658" align="aligncenter" width="1654"] Bàn đồ, sơ đồ tuyến điểm du lịch Việt Nam[/caption] Bản đồ du lịch Việt Nam (theo tỉnh/thành) Bản đồ du lịch Tây Bắc Vùng Tây Bắc  là vùng miền núi phía tây của  miền Bắc Việt Nam , có chung đường biên giới với Lào và Trung Quốc. Vùng này có khi được gọi là  Tây Bắc Bắc Bộ  và là một trong 3 tiểu vùng địa lý tự nhiên của Bắc Bộ Việt Nam (2 tiểu vùng kia là  Vùng Đông Bắc  và  Đồng bằng sông Hồng ). Về mặt hành chính, vùng Tây Bắ
Đọc thêm