Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Driving Distance
-- In case of commercial use please contact me
-- I don't mind earning money using my work if you share it with me
-- I am counting on your honesty
-- Function: public.ak_pgr_isochrone(double precision, double precision, double precision, double precision, integer, integer)
-- DROP FUNCTION public.ak_pgr_isochrone(double precision, double precision, double precision, double precision, integer, integer);
CREATE OR REPLACE FUNCTION public.ak_aba_pgr_isochrone(
start_x double precision,
start_y double precision,
driving_time double precision,
concave double precision,
epsg integer,
debug integer)
RETURNS geometry AS
$BODY$
DECLARE
start_time TIMESTAMP WITHOUT TIME ZONE;
start_x_4326 double precision;
start_y_4326 double precision;
start_edge record;
start_edge_i integer;
driving_time_h double precision;
res_geom geometry;
g_range geometry;
BEGIN
/************************************************************************************
*************************************************************************************
*** Funkcja rysująca izochrony (obszar dojazdu w zadanym czasie) ***
*** Funkcja szuka odcinków które w całości są dostępne w zadanym czasie ***
*** po czym znajduje wszystkie odcinki styczne z nimi i dodaje taki ich fragment ***
*** jaki można przejechać w czasie który pozostał. ***
*** Geomteria tworzona jest przez dodanie do każdego odcinka bufora o wartości ***
*** równej odległości możliwej do przejechania z prędkością 7km/h. ***
*** Na razie funkcja nie bierze pod uwagę przeszkód terenowych ponieważ ***
*** uwzględnienie ich znacząco wpływa na wydajność (wydłuża czas 10x) ***
*** ***
*** Parametry wejściowe: ***
*** - start_x - współrzędna x punktu początkowego ***
*** - start_y - współrzędna y punktu początkowego ***
*** - driving_time - czas graniczny w minutach ***
*** - concave double - współczynnik zassania (nieużywany) ***
*** - epsg - układ odniesienia współrzędnych (i zwracanej geometrii) ***
*** - debug - włączenie powoduje wyświetlenie dodatkowych komunikatów ***
*** ***
*** AK - luty 2016 ***
*************************************************************************************
************************************************************************************/
-- Przykładowe wywołanie funkcji dla KM PSP w BB
-- select * from aba_pgr_isochrone_new_test(19.0496201,49.8025666,59,0.7,4326,1);
if debug = 0 then
SET client_min_messages TO warning;
else
SET client_min_messages TO debug;
end if;
-- obsługa błędów parametrów wejściowych funkcji
if epsg = 4326 and (start_x < 14 or start_x > 24.2) then raise warning 'Długość geograficzna punktu startowego musi zawierać się w zakresie od 14 do 24.2 stopni'; end if;
if epsg = 4326 and (start_y < 49 or start_y > 56) then raise warning 'Szerokość geograficzna punktu startowego musi zawierać się w zakresie od 49 do 56 stopni'; end if;
if debug not in (0,1) then raise warning 'Parametr debug może przyjmować wartości 0 lub 1'; end if;
-- dopisać kontrolę dla innych epsg
--przeliczanie parametrów wejściowych na 4326
if epsg = 4326 then
start_x_4326 = start_x;
start_y_4326 = start_y;
if debug = 1 then raise notice 'Parametry wejściwe w EPSG:4326 - przeliczanie zbędne.';end if;
else
start_x_4326 = st_x(st_transform(ST_SetSRID(ST_MakePoint(start_x,start_y),epsg),4326));
start_y_4326 = st_y(st_transform(ST_SetSRID(ST_MakePoint(start_x,start_y),epsg),4326));
if debug = 1 then raise notice 'Parametry wejściwe po przeliczeniu: start(%,%)',start_x_4326, start_y_4326;end if;
end if;
-- recheck parametrów wejściowych po przeliczeniu
if start_x_4326 < 14 or start_x_4326 > 24.2 then
raise warning 'Długość geograficzna punktu startowego po przeliczeniu musi zawierać się w zakresie od 14 do 24.2 stopni';
return null;
end if;
if start_y_4326 < 49 or start_y_4326 > 56 then
raise warning 'Szerokość geograficzna punktu startowego po przeliczeniu musi zawierać się w zakresie od 49 do 56 stopni';
return null;
end if;
-- przeliczanie czasu jazdy
driving_time_h = driving_time/60;
if debug = 1 then
start_time := clock_timestamp();
raise notice 'Wczytywanie zakresu danych';
end if;
-- zakładam 120km/h i rysuję okrą gdzie można zalecieć z tą prędkością - wynikowe izochrony na pewno będą mniejsze
g_range = ST_Transform(ST_Buffer(ST_Transform(ST_SetSRID(ST_MakePoint(start_x_4326,start_y_4326),4326),2180),(2000*driving_time)),4326);
create temporary table test_edges on commit drop as
select id, osm_id, osm_name, osm_meta, osm_source_id, osm_target_id, clazz, flags, source, target, km, kmh, cost, reverse_cost, x1, y1, x2, y2, geom_way from pgr_poland where geom_way && g_range;
create index test_edges_geom_idx on test_edges using gist(geom_way);
if debug = 1 then
raise notice 'Zakres danych %, czas: %', (select count(1) from test_edges), extract(epoch from (clock_timestamp() - start_time))*1000;
end if;
if debug = 1 then
start_time := clock_timestamp();
raise notice 'Wyszukiwanie początku trasy';
end if;
-- szukamy początku trasy
SELECT source, ST_Distance(geom_way, ST_SetSRID(ST_MakePoint(start_x_4326,start_y_4326),4326)) AS dist
FROM test_edges WHERE geom_way && ST_SetSRID(ST_MakeBox2D(ST_MakePoint((start_x_4326-0.004),(start_y_4326-0.004)),ST_MakePoint((start_x_4326+0.002),(start_y_4326+0.002))), 4326)
ORDER BY dist LIMIT 1
INTO start_edge;
if not found then --wyjątek jak nie znajdzie
Raise warning 'Nie odnaleziono drogi w odległości 200m od punktu startowego';
end if;
start_edge_i = start_edge.source;
if debug = 1 then
raise notice 'Punkt poczatkowy %, czas: %', start_edge_i, extract(epoch from (clock_timestamp() - start_time))*1000;
end if;
if debug = 1 then
start_time := clock_timestamp();
raise notice 'Wyznaczenie czasów dostępności odcinków';
end if;
-- wyznaczam czas dostępności każdego odcinka w tabeli z pobranymi danymi
create temporary table test_edges_cost on commit drop as
SELECT row_number() over() as id, e.cost as cost_edge, r.agg_cost as cost_route, e.clazz, geom_way as geom, null::geometry as buffer
FROM test_edges e
JOIN (SELECT * FROM pgr_drivingdistance('select id, source, target, cost from test_edges', start_edge_i, 10, false))
AS r
ON e.source = r.node;
if debug = 1 then
raise notice 'Odcinki z obliczonym czasem: %, czas: %',(select count(1) from test_edges_cost),extract(epoch from (clock_timestamp() - start_time))*1000;
end if;
if debug = 1 then
start_time := clock_timestamp();
raise notice 'Wznaczenie odcinków w całości w izochronie';
end if;
-- wyciągam z tabeli odcinki, które w całości mieszczą się w założonym czasie
create temporary table test_edges_final on commit drop as
select * from test_edges_cost c where c.cost_route <= driving_time_h;
create index test_edges_final_geom_idx on test_edges_final using gist(geom);
create index test_edges_final_clazz_idx on test_edges_final using btree(clazz);
create index test_edges_final_id_idx on test_edges_final using btree(id);
create index test_edges_final_cost_route_idx on test_edges_final using btree(cost_route);
if debug = 1 then
raise notice 'Całych odcinków: %, czas: %',(select count(1) from test_edges_final),extract(epoch from (clock_timestamp() - start_time))*1000;
end if;
if debug = 1 then
start_time := clock_timestamp();
raise notice 'Wznaczenie odcinków częściowo w izochronie';
end if;
-- wyciągam z tabeli odcinki, które przecinają się z tymi wybranymi do izochrony
create temporary table test_edges_border on commit drop as
select c.*, ef.id as main_id, (driving_time_h - ef.cost_route) as left_cost, null::geometry as blade from test_edges_cost c
join test_edges_final ef on ST_Intersects(ef.geom, c.geom) and c.cost_route > driving_time_h
where (driving_time_h - ef.cost_route) < c.cost_edge;
create index test_edges_border_main_id_idx on test_edges_border using btree(main_id);
if debug = 1 then
raise notice 'Częściowych odcinków: %, czas: %',(select count(1) from test_edges_border),extract(epoch from (clock_timestamp() - start_time))*1000;
end if;
if debug = 1 then
start_time := clock_timestamp();
raise notice 'Processing odcinków częściowych';
end if;
-- wyszukuje odcinki biegnące 'do izochrony' i odwracam ich kierunek
update test_edges_border ef
set geom = ST_Reverse(ef.geom)
where ST_Intersects(ST_LineInterpolatePoint(geom,0), (select geom from test_edges_final c where c.id=ef.main_id))=false;
if debug = 1 then
raise notice 'Odwracanie odcinków: %',extract(epoch from (clock_timestamp() - start_time))*1000;
start_time := clock_timestamp();
end if;
-- usuwam duble odcinków stycznych do izochrony
delete from test_edges_border b1 where exists (select 1 from test_edges_border b2 where b1.id = b2.id and b1.left_cost > b2.left_cost);
if debug = 1 then
raise notice 'Usuwanie dubli: %',extract(epoch from (clock_timestamp() - start_time))*1000;
start_time := clock_timestamp();
end if;
--wyznaczam punkty cięcia odcinków dodatkowych
update test_edges_border set blade = ST_LineInterpolatePoint(geom,(left_cost/cost_edge));
if debug = 1 then
raise notice 'Wyznaczenie punktów cięcia: %',extract(epoch from (clock_timestamp() - start_time))*1000;
start_time := clock_timestamp();
end if;
-- przycinam odcinki dodatkowe usuwając ich nadmiarowe części
create temporary table test_edges_border_cutted on commit drop as
select aa.* from
(select id, cost_edge, cost_route, buffer, main_id, left_cost, (ST_Dump(ST_Difference(geom, ST_Transform(ST_Buffer(ST_Transform(blade,2180),1),4326)))).geom as geom
from test_edges_border) as aa
join test_edges_final f on f.id = aa.main_id and ST_Intersects(aa.geom, f.geom) ;
if debug = 1 then
raise notice 'Cięcie odcinków częściowych z rozbiciem: %',extract(epoch from (clock_timestamp() - start_time))*1000;
start_time := clock_timestamp();
end if;
-- zakładam bufor dla odcinków dodatkowych (stały - 10m)
update test_edges_border_cutted set buffer = ST_Transform(ST_Buffer(ST_Transform(geom,2180),10),4326);
if debug = 1 then
raise notice 'Wyznaczenie bufora odcinków częściowych: %, czas: %',(select count(1) from test_edges_border_cutted), extract(epoch from (clock_timestamp() - start_time))*1000;
start_time := clock_timestamp();
end if;
-- wyznaczam bufor dla odcinków pełnych (z wyjątkiem A i S) o wartości pozostałego czasu dla odcinka
update test_edges_final set buffer = ST_Transform(ST_Buffer(ST_Transform(geom,2180),((driving_time_h-cost_route)*60*100)+10),4326) where clazz > 14;
if debug = 1 then
raise notice 'Wyznaczenie bufora odcinków pełnych: %',extract(epoch from (clock_timestamp() - start_time))*1000;
start_time := clock_timestamp();
end if;
-- wyznaczam bufor dla odcinków S i A wartości 20m (te drogi są ogrodzone)
update test_edges_final set buffer = ST_Transform(ST_Buffer(ST_Transform(geom,2180),20),4326) where clazz <=14;
if debug = 1 then
raise notice 'Wyznaczenie bufora Autostrad i ekspresów: %',extract(epoch from (clock_timestamp() - start_time))*1000;
start_time := clock_timestamp();
end if;
-- łącze otrzymane bufory w jeden
res_geom =
(
select ST_Union(buffer)
from (
select buffer from test_edges_border_cutted
where ST_IsValid(buffer)
union all
select buffer from test_edges_final
where ST_IsValid(buffer)
) as aa
);
if debug = 1 then
raise notice 'Połączenie buforów: %', extract(epoch from (clock_timestamp() - start_time))*1000;
end if;
-- zwracam geometrię z konwersją do żądanego układu współrzędnych
return st_transform(res_geom, epsg);
END;
$BODY$
LANGUAGE plpgsql VOLATILE
COST 100;
ALTER FUNCTION public.ak_pgr_isochrone(double precision, double precision, double precision, double precision, integer, integer)
OWNER TO osm;
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.