update
denver_buildings
set
geom_z = st_force3d(
-- First value is the building geometry in EPSG 4326
st_transform(geom, 4326),
-- This is our Z or height value, the building height
-- plus the ground elevation
bldg_heigh + ground_ele
)
-- Find our first building
with a as (
select
geom_z,
gid
from
denver_buildings
limit
1 offset 100
),
-- Find a building with its ID and GeometryZ
-- within 2 kilometers
b as (
select
geom_z,
gid
from
denver_buildings
where
st_dwithin(
st_transform(geom, 3857),
(
select
st_transform(geom, 3857)
from
a
),
2000
)
limit
1
),
-- Use UNION to create a single table with matching columns
c as (
select
*
from
a
union
select
*
from
b
),
-- Store the geometries and and IDs in arrays
bldgs as (
select
array_agg(st_centroid(geom_z)) as geom,
array_agg(gid) as id
from
c
),
-- This query finds all the buildings within 3 kilometers of each building
denver as (
select
st_transform(geom, 3857) as geom,
geom_z,
gid
from
denver_buildings
where
st_dwithin(
-- We union the two points and turn them back into 2D
-- Geometries so we can run the query to find all
-- buildings with 3 kilometers one time
st_union(
(
select
st_force2d(st_transform(geom [1], 3857))
from
bldgs
),
(
select
st_force2d(st_transform(geom [2], 3857))
from
bldgs
)
),
st_transform(geom, 3857),
3000
)
)
select
d.gid,
st_transform(d.geom, 4326) as geom
from
denver d
-- Now we can use our ST_3DIntersects funciton to find
-- all the buildings crossing our path
join bldgs on st_3dintersects(
-- The path can be built using the same ST_MakeLine
-- function we have been using
st_makeline(bldgs.geom [1], bldgs.geom [2]),
d.geom_z
)