Using PostGIS for isovists calculation
First of all, what is an isovist
? It’s simply the same as a viewshed
… and what the heck is a viewshed
? Let’s go to the Wikipedia:
A single isovist is the volume of space visible from a given point in space,
together with a specification of the location of that point
The decision of using the word isovist
or viewshed
is only related to the knowledge field you’re working on:
 Isovist: architecture, space syntax, landscaping, urbanism, geomarketing
 viewshed: geography, infrastructure location
In this article, I’m going to talk about the 2D version of the problem, but could be expanded to cover the 3D case.
Let’s define a real use case as a base for later explanations. Let’s say we’re and ad company and want to place billboards in the facades of buildings in the center of Madrid, targeting the people waiting for the bus at the defined stops. Our target would be to find the facades with the highest eyeballs
audience. The result should be like this
Cool, isn’t it? And how is it done? The answer is quite simple
Yes, magic.
And PostGIS.
Ok, follow me in this journey along this raytracinglike approach using a somehow tricky PL/pgSQL function
Let’s say, that the point I’m going to study (center
) has a heading of 270º and a fov of 90º, and the horizon is located at 150m. The arbitrary point for the following images is defined as:
SELECT ST_SetSRID(ST_MakePoint( 3.705856,40.420371),4326) as center;
First of all, let’s think on the macro structure of the problem. We will need a point geometry for any input point and a set of buildings
polygons. But… we will need some extra parameters too, like:
radius
: Horizon distance, so we don’t overcharge the calculation with useless inforays
: number of angular steps, that will define the resolutionheading
: northing in degrees, to set the direction theeyeballs
are looking atfov
: field of view in degrees, centered in the heading direction
So the signature of the function would be like:
CREATE OR REPLACE FUNCTION ISOVIST(
IN center geometry,
IN polygons geometry[],
IN radius numeric DEFAULT 150,  horizon distance in meters
IN rays integer DEFAULT 36,  number of rays
IN heading integer DEFAULT 999,  take heading (degrees, 0360) into account
IN fov integer DEFAULT 360  field of view (degrees)
)
Having a default value for fov
of 360º, so the heading is not needed in this very case.
Once there, let’s declare some support variables
arc numeric;  angular resolution
angle_0 numeric;  starting angle for fov != 360º
And the output
geomout geometry;  resulting geometry of the isovist
And calc the values for the support variables:
 resolution in degrees
arc := fov::numeric / rays::numeric;
 fov range start value
IF fov = 360 THEN
angle_0 := 0;
ELSE
 centered in heading angle
angle_0 := heading  0.5 * fov;
END IF;
Now, we’re ready to start querying the data. We’re building it with incremental subqueries using WITH
clause, so each block can be easily explained:

Unnest the polygons array to have a table object:
buildings_0 AS( SELECT t.geom FROM unnest(polygons) as t(geom) ),

Filter out the polygons further than the horizon, so we limit the next calculi to the polygons of interest only
buildings_crop AS( SELECT geom FROM buildings_0 WHERE ST_DWithin( center::geography, geom::geography, radius ) ),

Let’s add the horizon as
the final polygon
that stops any ray traced fromcenter
buildings AS( SELECT geom FROM buildings_crop UNION ALL SELECT ST_buffer(center::geography, radius)::geometry as geom ),

Now, we need to create #
rays
rays, fromcenter
, eacharc
degrees starting atangle_0
. To do so, we’re gonna use the PostGIS function ST_Project that generates a 2nd point from a starting one, a direction and a distance. So my nth ray will look likeST_SetSRID( ST_MakeLine( center, ST_Project( center::geography, radius + 1, radians(angle_0 + n::numeric * arc) )::geometry ), 4326) AS geom
We use
radius+1
as length to be sure that the ray intersects the horizon at least.Using generate_series we can create all the rays at once:
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) ),

Now, getting the intersections between the rays and the polygons:
intersections AS( SELECT r.id, (ST_Dump( ST_Intersection(  to avoid intersections with polygon holes ST_Boundary(b.geom), r.geom ) )).geom AS point FROM rays r LEFT JOIN buildings b ON ST_Intersects(b.geom,r.geom) ),
So we translated our rays spectrum into a cloud of points of the intersections between all the rays and all the (external boundaries of the) polygons within range

Now, we need to rank the intersections in terms of distance to
center
…intersections_distances AS( SELECT id, point as geom, row_number() over(partition by id order by center <> point) as ranking FROM intersections ),

…in order to keep only the closest intersection per polygon
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 ),
We are adding the center point twice (before and after all the intersections) in case the
fov
is less than 360º… 
… so we can close the line and build a polygon (the isovist!) with it
isovist_0 AS( SELECT ST_MakePolygon(ST_MakeLine(geom)) as geom FROM intersection_closest ),

Let’s get the
polygons
that actually intersects my isovistisovist_buildings AS( SELECT  avoid geometry collections ST_CollectionExtract(ST_union(b.geom),3) as geom FROM isovist_0 i, buildings_crop b WHERE ST_Intersects(b.geom,i.geom) )

Because of the
arc
resolution, we need to crop the resulting isovist with thepolygons
above to have the best level of detail. Because of this PostGIS ticket we need to sanitize the results for the cases there’s no building intersecting my isovist (all of them are further than the horizon), and it’s a plain circular arc.
SELECT
coalesce(ST_Difference(i.geom, b.geom), i.geom) into geomout
FROM
isovist_0 i,
isovist_buildings b;
So, finally, the full function should look like:
CREATE OR REPLACE FUNCTION ISOVIST(
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;
angle_0 numeric;
geomout geometry;
BEGIN
arc := fov::numerics / 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;
Magic!