Driving Distance

gistfile1.txt

-- 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;