Using PostGIS for isovists calculation

6 minute read

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 ray-tracing-like 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 info
rays: number of angular steps, that will define the resolution
heading: northing in degrees, to set the direction the eyeballs are looking at
fov: 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, 0-360) 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 from center

    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, from center, each arc degrees starting at angle_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 n-th ray will look like

ST_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 isovist

    isovist_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 the polygons 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!

Share on

Twitter Facebook Google+ LinkedIn

Abel Vázquez

Using PostGIS for isovists calculation

Share on

You May Also Enjoy

Fast contour lines using plain PostGIS

Some little known PostgreSQL/PostGIS tips & tricks

QuadGrid: fast spatial data generalization for the best LoD

Dymanically find loops in a tracks dataset