· 6 years ago · Jun 19, 2019, 06:02 AM
1-- Function: public.update_visibility_graph(double precision, double precision)
2
3-- DROP FUNCTION public.create_visibility_graph(text, text, double precision, double precision);
4
5CREATE OR REPLACE FUNCTION public.create_visibility_graph(
6
7IN barriers text,
8
9IN sites text,
10
11IN max_distance double precision DEFAULT 40000,
12
13IN buffer_distance double precision DEFAULT 100,
14
15OUT result text)
16
17RETURNS text AS
18
19$BODY$
20
21DECLARE
22
23querystring text;
24
25BEGIN
26
27drop table if exists public.visibility_graph;
28
29drop table if exists public.allpts;
30
31-- create a table for the points
32
33create temporary table allpts (
34
35gid serial primary key,
36
37id integer,
38
39the_geom geometry(Point,3857) not null
40
41);
42
43-- populate the temp table with convex coastal points
44
45querystring = $$
46
47insert into allpts (the_geom)
48
49select
50
51distinct the_geom
52
53from (
54
55select the_geom,
56
57(angle > 180) or (angle between -180 and 0) as concave
58
59from (
60
61select pt as the_geom,
62
63degrees(atan2((st_y(ptn) - st_y(pt)), (st_x(ptn) - st_x(pt)))) - degrees(atan2((st_y(ptp) - st_y(pt)), (st_x(ptp) - st_x(pt)))) as angle
64
65from (
66
67select
68
69pt, feat, ptnum,
70
71lag(pt,1) over (partition by feat) as ptp, -- the previous point, or the last if null lag(pt,-1,pt2) over (partition by feat) as ptn -- the next point or null if it is the last
72
73from (
74
75select
76
77(dmp).path[1] as feat,
78
79(st_dumppoints((dmp).geom)).geom as pt,
80
81(st_dumppoints((dmp).geom)).path[2] as ptnum, st_pointn(ST_ExteriorRing((dmp).geom),2) as pt2
82
83from (
84
85select
86
87st_dump(
88
89st_simplifypreservetopology(
90
91st_buffer(
92
93st_intersection(barriers.geom, sitebuffer.geom), $$||buffer_distance||$$, 1) ,$$||buffer_distance||$$)
94
95) as dmp from (
96
97-- get the barrier polygon layer
98
99select st_transform(st_union(geom),3857) as geom
100
101from $$||barriers||$$) barriers
102
103join (
104
105-- create a buffer of the sites
106
107select st_union(st_buffer(st_transform(geom,3857), $$||max_distance||$$)) as geom from $$||sites||$$) sitebuffer on true
108
109) d
110
111) foo
112
113) bar where ptp is not null -- exclude the first point of each polygon (same as last) ) foobar
114
115join m50nad17n b on not st_intersects(b.geom, the_geom)
116
117) foofoobar where not concave$$;
118
119execute querystring;
120
121-- add in the site locations
122
123querystring := $$insert into allpts (id, the_geom)
124
125select s.gid, st_transform(s.geom, 3857)
126
127from $$||sites||$$ s$$;
128
129execute querystring;
130
131-- join all nearby point pairs that have a clear line of sight
132
133querystring:=$$create table visibility_graph as
134
135select
136
137row_number() over () as id,
138
139p1.gid as source,
140
141p2.gid as target,
142
143p1.id as sid, p2.id as tid,
144
145st_makeline(p1.the_geom, p2.the_geom)::geometry(LINESTRING,3857) as the_geom
146
147from allpts p1
148
149join allpts p2 on p1.gid < p2.gid
150
151and st_distance(p1.the_geom, p2.the_geom) < max_distance
152
153join $$||barriers||$$ l2 on true -- obstacle layer
154
155where not st_intersects(l2.geom, st_makeline(p1.the_geom, p2.the_geom))$$;
156
157execute querystring;
158
159-- done with the temp table drop table allpts;
160
161-- build routing table
162
163select into result pgr_createTopology('visibility_graph',0.0001, clean:=true);
164
165END
166
167$BODY$
168
169LANGUAGE
170
171plpgsql VOLATILE
172
173COST 100;
174
175
176
177
178
179SET search_path TO public;
180
181create_visibility_graph(m50nad17n, sites2, 40000, 100);