set search_path to d16, osm, ruian, public;

----------------------------------------------------------------------------------------------
-- silnice

-- vyber silnic
create table roads as
select	highway as categor,
		oneway,
		st_transform(geom, 5514) as geom
from	czech_line
where	highway in
		(	'motorway', 'trunk', 'primary', 'secondary', 'tertiary',
			'motorway_link', 'trunk_link', 'primary_link', 'secondary_link', 'tertiary_link');
select updateGeometrySRID('roads', 'geom', 5514);

-- unikatni typy geometrie
select	distinct st_geometryType(geom)
from	roads;

-- primarni klic a prostorovy index
alter table roads add column gid serial;
alter table roads add primary key (gid);
create index on roads using gist (geom);

-- sloupce pro pgRouting
alter table roads add column source integer;
alter table roads add column target integer;
create index on roads(source);
create index on roads(target);

-- sloupec s delkou
alter table roads add column length double precision;
update roads set length = st_length(geom);

-- sloupec s rychlosti
alter table roads add column speed integer;
update roads
set speed = case categor
				when 'motorway' then 110
				when 'motorway_link' then 85
				when 'trunk' then 105
				when 'trunk_link' then 80
				when 'primary' then 65
				when 'primary_link' then 50
				when 'secondary' then 50
				when 'secondary_link' then 40
				when 'tertiary' then 40
				when 'tertiary_link' then 30
			end;

-- sloupec s casem
alter table roads add column time double precision;
update roads set time = length/(speed/3.6);

-- sloupce pro jednosmerne komunikace
alter table roads add column this_way double precision;
update roads
set this_way =	case oneway
					when '-1' then 1000000000
					when 'reverse' then 1000000000
					else 0
				end;

alter table roads add column oppo_way double precision;
update roads
set oppo_way =	case oneway
					when 'yes' then 1000000000
					when 'true' then 1000000000
					when '1' then 1000000000
					else 0
				end;

-- vytvoreni topologie
select pgr_createTopology('roads', 0.001, 'geom', 'gid');

-- analyza vrstvy
select pgr_analyzeGraph('roads', 0.001, 'geom', 'gid', 'source', 'target');

-- oprava topologie
select pgr_nodeNetwork('roads', 0.001, 'gid', 'geom', 'noded');

-- vytvoreni topologie
select pgr_createTopology('roads_noded', 0.001, 'geom', 'id');

-- analyza vrstvy
select pgr_analyzeGraph('roads_noded', 0.001, 'geom', 'id', 'source', 'target');

-- oprava topologie 2
select pgr_nodeNetwork('roads_noded', 0.001, 'id', 'geom', '2');

-- vytvoreni topologie 2
select pgr_createTopology('roads_noded_2', 0.001, 'geom', 'id');

-- analyza vrstvy 2
select pgr_analyzeGraph('roads_noded_2', 0.001, 'geom', 'id', 'source', 'target');

-- vytvoreni finalni vrstvy
create table roads_f as
select		rn2.id as gid, r.categor as categor, r.oneway as oneway,
			rn2.geom as geom
from		roads_noded_2 rn2
left join	roads_noded rn on rn2.old_id = rn.id
left join	roads r on rn.old_id = r.gid;

-- primarni klic a prostorovy index
alter table roads_f add primary key (gid);
create index on roads_f using gist (geom);

-- sloupce pro pgRouting
alter table roads_f add column source integer;
alter table roads_f add column target integer;
create index on roads_f(source);
create index on roads_f(target);

-- sloupec s delkou
alter table roads_f add column length double precision;
update roads_f set length = st_length(geom);

-- sloupec s rychlosti
alter table roads_f add column speed integer;
update roads_f
set speed = case categor
				when 'motorway' then 110
				when 'motorway_link' then 85
				when 'trunk' then 105
				when 'trunk_link' then 80
				when 'primary' then 65
				when 'primary_link' then 50
				when 'secondary' then 50
				when 'secondary_link' then 40
				when 'tertiary' then 40
				when 'tertiary_link' then 30
			end;

-- sloupec s casem
alter table roads_f add column time double precision;
update roads_f set time = length/(speed/3.6);

-- sloupce pro jednosmerne komunikace
alter table roads_f add column this_way double precision;
update roads_f
set this_way =	case oneway
					when '-1' then 1000000000
					when 'reverse' then 1000000000
					else 0
				end;

alter table roads_f add column oppo_way double precision;
update roads_f
set oppo_way =	case oneway
					when 'yes' then 1000000000
					when 'true' then 1000000000
					when '1' then 1000000000
					else 0
				end;

-- vytvoreni topologie
select pgr_createTopology('roads_f', 0.001, 'geom', 'gid');

-- kontrola geometrie
select	gid
from	roads_f
where	st_isValid(geom) != true;

----------------------------------------------------------------------------------------------
-- vesnice
create table villages as
select	osm_id as gid, name, st_transform(geom, 5514) as geom
from	czech_point where osm_id in (359094467, 1599445934);

select updateGeometrySRID('villages', 'geom', 5514);
alter table villages add primary key (gid);
create index on villages using gist (geom);

-- kontrola geometrie
select	gid
from	villages
where	st_isValid(geom) != true;

----------------------------------------------------------------------------------------------
-- reky
create table rivers as
select		name, st_transform(st_collect(geom), 5514) as geom
from		czech_line
where 		waterway = 'river'
			and
			name is not null
group by	name;

select updateGeometrySRID('rivers', 'geom', 5514);
alter table rivers add column gid serial;
alter table rivers add primary key (gid);
create index on rivers using gist (geom);

-- kontrola geometrie
select	gid
from	rivers
where	st_isValid(geom) != true;

----------------------------------------------------------------------------------------------
-- letiste
create table airports as (
select		name, st_transform(st_centroid(st_collect(geom)), 5514) as geom
from		(	select		name, st_centroid(geom) as geom
				from		czech_polygon
				where		aeroway = 'aerodrome'
							and
							name is not null
				union
				select		name, geom
				from		czech_point
				where		aeroway = 'aerodrome'
							and
							name is not null) a
group by	name);

select updateGeometrySRID('airports', 'geom', 5514);
alter table airports add column gid serial;
alter table airports add primary key (gid);
create index on airports using gist (geom);

-- kontrola geometrie
select	gid
from	airports
where	st_isValid(geom) != true;

----------------------------------------------------------------------------------------------
-- priklad 1

-- nejrychlejsi Dijkstra
create table path_d_time as
select		seq, path_seq, node, edge, cost, agg_cost, b.length, b.geom as geom
from	 	pgr_dijkstra(
			'select gid as id, source, target, (time+this_way) as cost, (time+oppo_way) as reverse_cost from roads_f',
			(	select	id
				from 	roads_f_vertices_pgr r,
						(select geom from villages where name = 'Onen Svět') os
				order by st_distance(r.the_geom, os.geom) limit 1),
			(	select	id
				from	roads_f_vertices_pgr as r,
						(select geom from villages where name = 'Peklo') p
				order by st_distance(r.the_geom, p.geom) limit 1)) a
left join	roads_f b on a.edge = b.gid;

-- delka
select sum(length)/1000 as length from path_d_time;
-- 307.8 km

--cas
select agg_cost/60 as time from path_d_time where edge = -1;
-- 211 minut

-- nejkratsi Dijkstra
create table path_d_len as
select		seq, path_seq, node, edge, cost, agg_cost, b.time, b.geom as geom
from	 	pgr_dijkstra(
			'select gid as id, source, target, (length+this_way) as cost, (length+oppo_way) as reverse_cost from roads_f',
			(	select	id
				from 	roads_f_vertices_pgr r,
						(select geom from villages where name = 'Onen Svět') os
				order by st_distance(r.the_geom, os.geom) limit 1),
			(	select	id
				from	roads_f_vertices_pgr as r,
						(select geom from villages where name = 'Peklo') p
				order by st_distance(r.the_geom, p.geom) limit 1)) a
left join	roads_f b on a.edge = b.gid;

-- delka
select agg_cost/1000 as length from path_d_len where edge = -1;
-- 294.8 km

--cas
select sum(time)/60 as time from path_d_len;
-- 340 min

-- sloupce pro A*
alter table roads_f add column x1 double precision;
update roads_f set x1 = st_x(st_startpoint(geom));

alter table roads_f add column y1 double precision;
update roads_f set y1 = st_y(st_startpoint(geom));

alter table roads_f add column x2 double precision;
update roads_f set x2 = st_x(st_endpoint(geom));

alter table roads_f add column y2 double precision;
update roads_f set y2 = st_y(st_endpoint(geom));

----------------------------------------------------------------------------------------------
-- priklad 2

-- nejrychlejsi A*
create table path_a_len as
select		seq, id1 as node, id2 as edge, cost, b.time, b.geom as geom
from	 	pgr_astar(
			'select gid as id, source, target, (length+this_way) as cost, x1, y1, x2, y2, (length+oppo_way) as reverse_cost from roads_f',
			cast((	select	id from roads_f_vertices_pgr r,
					(select geom from villages where name = 'Onen Svět') os
					order by st_distance(r.the_geom, os.geom) limit 1) as int),
			cast((	select id from roads_f_vertices_pgr r,
					(select geom from villages where name = 'Peklo') p
					order by st_distance(r.the_geom, p.geom) limit 1) as int),
			true, true) a
left join	roads_f b on a.id2 = b.gid;

-- delka
select sum(cost)/1000 as length from path_a_len;
-- 294.8 km

-- cas
select sum(time)/60 as time from path_a_len;
-- 340 min

----------------------------------------------------------------------------------------------
-- priklad 3

-- projete obce
with path_d_len_o as (
	select	distinct 	o.kod, o.geom
	from	path_d_len p
	join	obce o on	p.geom && o.geom
						and
						st_intersects(p.geom, o.geom))
select	count(*)
from	path_d_len_o;
-- 84 obci

----------------------------------------------------------------------------------------------
-- priklad 4

-- obec s nejvetsi najetou vzdalenosti
with path_d_len_o as (
	select	distinct o.nazev, o.geom
	from	path_d_len p
	join	obce o on	p.geom && o.geom
						and
						st_intersects(p.geom, o.geom))
select		o.nazev,
			(st_length(st_intersection(o.geom, u.geom)))/1000 as length
from		path_d_len_o o
join (		select	st_collect(geom) as geom
			from	path_d_len) u on o.geom && u.geom
order by	length desc limit 1;
-- Pardubice, 14.3 km

----------------------------------------------------------------------------------------------
-- priklad 5

-- prekrocene reky
select	ri.name, count(*)
from	path_d_len p
join	rivers ri on	p.geom && ri.geom
						and
						st_crosses(p.geom, ri.geom)
group by	ri.name
having	count(*) > 1;
-- Otava, 3
-- Zdobnice, 2

----------------------------------------------------------------------------------------------
-- priklad 6

-- doba cesty pri opousteni Zapadoceskeho kraje
with path_d_len_i as (
	select		p.seq, p.cost as length, r.speed, st_intersection(p.geom, k.geom), k.kod as code
	from		path_d_len p
	left join (	select	kod, geom
				from	kraje
				where	nazev = 'Západočeský') k
				on	p.geom && k.geom
					and
					st_intersects(p.geom, k.geom)
	left join	roads_f r on p.edge = r.gid)
select	sum(length/(speed/3.6))/60 as time
from	path_d_len_i
where	seq < (	select	min(seq)
				from	path_d_len_i
				where	code is null);
-- 61 min

----------------------------------------------------------------------------------------------
-- priklad 7

-- pocet letist ve vzdalenosti 2 km od cesty
with	path_d_len_b as (
			select	st_buffer(geom, 2000) as geom
			from (	select	st_collect(geom) as geom
					from	path_d_len) u),
		path_d_time_b as (
			select	st_buffer(geom, 2000) as geom
			from (	select	st_collect(geom) as geom
					from	path_d_time) u)
select	'shortest' as path_type, count(*) as airports
from	airports a
join	path_d_len_b b on	a.geom && b.geom
							and
							st_within(a.geom, b.geom)
union
select	'fastest' as path_type, count(*) as airports
from	airports a
join	path_d_time_b b on	a.geom && b.geom
							and
							st_within(a.geom, b.geom);
-- fastest, 3
-- shortest, 8

