-- Function: public.update_visibility_graph(double precision, double precision)

-- DROP FUNCTION public.create_visibility_graph(text, text, double precision, double precision);

CREATE OR REPLACE FUNCTION public.create_visibility_graph(

IN barriers text,

IN sites text,

IN max_distance double precision DEFAULT 40000,

IN buffer_distance double precision DEFAULT 100,

OUT result text)

RETURNS text AS

$BODY$

DECLARE

querystring text;

BEGIN

drop table if exists public.visibility_graph;

drop table if exists public.allpts;

-- create a table for the points

create temporary table allpts (

gid serial primary key,

id integer,

the_geom geometry(Point,3857) not null

);

-- populate the temp table with convex coastal points

querystring = $$

insert into allpts (the_geom)

select

distinct the_geom

from (

select the_geom,

(angle > 180) or (angle between -180 and 0) as concave

from (

select pt as the_geom,

degrees(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

from (

select

pt, feat, ptnum,

lag(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

from (

select

(dmp).path[1] as feat,

(st_dumppoints((dmp).geom)).geom as pt,

(st_dumppoints((dmp).geom)).path[2] as ptnum, st_pointn(ST_ExteriorRing((dmp).geom),2) as pt2

from (

select

st_dump(

st_simplifypreservetopology(

st_buffer(

st_intersection(barriers.geom, sitebuffer.geom), $$||buffer_distance||$$, 1) ,$$||buffer_distance||$$)

) as dmp from (

-- get the barrier polygon layer

select st_transform(st_union(geom),3857) as geom

from $$||barriers||$$) barriers

join (

-- create a buffer of the sites

select st_union(st_buffer(st_transform(geom,3857), $$||max_distance||$$)) as geom from $$||sites||$$) sitebuffer on true

) d

) foo

) bar where ptp is not null -- exclude the first point of each polygon (same as last) ) foobar

join m50nad17n b on not st_intersects(b.geom, the_geom)

) foofoobar where not concave$$;

execute querystring;

-- add in the site locations

querystring := $$insert into allpts (id, the_geom)

select s.gid, st_transform(s.geom, 3857)

from $$||sites||$$ s$$;

execute querystring;

-- join all nearby point pairs that have a clear line of sight

querystring:=$$create table visibility_graph as

select

row_number() over () as id,

p1.gid as source,

p2.gid as target,

p1.id as sid, p2.id as tid,

st_makeline(p1.the_geom, p2.the_geom)::geometry(LINESTRING,3857) as the_geom

from allpts p1

join allpts p2 on p1.gid < p2.gid

and st_distance(p1.the_geom, p2.the_geom) < max_distance

join $$||barriers||$$ l2 on true -- obstacle layer

where not st_intersects(l2.geom, st_makeline(p1.the_geom, p2.the_geom))$$;

execute querystring;

-- done with the temp table drop table allpts;

-- build routing table

select into result pgr_createTopology('visibility_graph',0.0001, clean:=true);

END

$BODY$

LANGUAGE

plpgsql VOLATILE

COST 100;





SET search_path TO public;

create_visibility_graph(m50nad17n, sites2, 40000, 100);