with buildings_0 as(
select
t.geom
from
-- Find all the buildingswithin 2 kilometers of Times Square where
-- geometries are stored in an array
unnest(
(
select
array_agg(geom)
from
nyc_building_footprints
where
st_dwithin(
geom :: geography,
st_setsrid(st_makepoint(-73.985136, 40.758786), 4326) :: geography,
2000
)
)
) as t(geom)
),
buildings_crop as(
-- Now only the buildings within 630 meters of Times Square
select
geom
from
buildings_0
where
st_dwithin(
st_setsrid(st_makepoint(-73.985136, 40.758786), 4326) :: geography,
geom :: geography,
630
)
),
buildings as(
-- Union these two tables together
select
geom
from
buildings_crop
union
all
select
st_buffer(
st_setsrid(st_makepoint(-73.985136, 40.758786), 4326) :: geography,
630
) :: geometry as geom
)
select
*
from
buildings
create
or replace function isovist(
-- Sets up the arguments for our function
in center geometry,
in polygons geometry [],
in radius numeric default 150,
in rays integer default 36,
in heading integer default -999,
in fov integer default 360
) returns geometry as $$ declare arc numeric;
-- Creates static variables for angle and geometry
angle_0 numeric;
geomout geometry;
-- Calculates the arc distance using the field of view (fov)
-- degrees and the number of rays
begin arc := fov :: numeric / rays :: numeric;
if fov = 360 then angle_0 := 0;
else angle_0 := heading - 0.5 * fov;
end if;
with buildings_0 as(
select
t.geom
from
unnest(polygons) as t(geom)
),
buildings_crop as(
select
geom
from
buildings_0
where
st_dwithin(center :: geography, geom :: geography, radius)
),
buildings as(
select
geom
from
buildings_crop
union
all
select
st_buffer(center :: geography, radius) :: geometry as geom
),
rays as(
select
t.n as id,
st_setsrid(
st_makeline(
center,
st_project(
center :: geography,
radius + 1,
radians(angle_0 + t.n :: numeric * arc)
) :: geometry
),
4326
) as geom
from
generate_series(0, rays) as t(n)
),
intersections as(
select
r.id,
(
st_dump(st_intersection(st_boundary(b.geom), r.geom))
).geom as point
from
rays r
left join buildings b on st_intersects(b.geom, r.geom)
),
intersections_distances as(
select
id,
point as geom,
row_number() over(
partition by id
order by
center <-> point
) as ranking
from
intersections
),
intersection_closest as(
select
-1 as id,
case
when fov = 360 then null :: geometry
else center
end as geom
union
all (
select
id,
geom
from
intersections_distances
where
ranking = 1
order by
id
)
union
all
select
999999 as id,
case
when fov = 360 then null :: geometry
else center
end as geom
),
isovist_0 as(
select
st_makepolygon(st_makeline(geom)) as geom
from
intersection_closest
),
isovist_buildings as(
select
st_collectionextract(st_union(b.geom), 3) as geom
from
isovist_0 i,
buildings_crop b
where
st_intersects(b.geom, i.geom)
)
select
coalesce(st_difference(i.geom, b.geom), i.geom) into geomout
from
isovist_0 i,
isovist_buildings b;
return geomout;
end;
$$ language plpgsql immutable;
select
*
from
isovist(
-- This is our center point in Times Square
st_setsrid(st_makepoint(-73.985136, 40.758786), 4326),
(
-- Selects all the geometries
select
array_agg(geom)
from
nyc_building_footprints
where
st_dwithin(
geom :: geography,
st_setsrid(st_makepoint(-73.985136, 40.758786), 4326) :: geography,
2000
)
),
300,
36
)