Mapping where the data is not

We discuss how to show areas on a web map where there is no data.

This technical item came up in a project where we needed to show three types of areas: areas that had a high flood risk, areas that had a low flood risk, and areas that had not been mapped for flood risk.

In this case we, we are using vector tiles drawn from a PostGIS database. We don’t want to constantly compute the complement of the data because that means unnecessary and prone to errors. Also, for many maps, the complement would be the rest of the world.

We’ll consider a simpler example for this post. We consider a simple PostGIS table containing the states of India. We want to create a map that displays the name of each Indian state when hovered over and displays “Out” when the mouse is outside of India.

We write the code in Clockwork Micro’s db2vector tool (https://www.clockworkmicro.com/services/map-maker), which is a UI tool to convert a geospatial sql query to vector tiles.

Let’s jump to what makes this query uniques. When a tile is fully inside India, each state polygon will have the name of the state as a feature and we can display these on hover. When a tile is fully outside of India, no data at all will be received from the database, yet we want to display the name “Out”. And when a tile partially intersects India (but not completely), the polygons within the intersection need to display the corresponding name field, while the complement of the Indian states within the tile needs to be assigned the name “Out”.

We have to take a peculiarity of ST_Difference into account: ST_Difference of a geometry and an empty statement is an empty statement rather than the geometry. And so we first have to check for intersection.

The segment of the code unique to our query is:

CASE
      WHEN ST_Intersects(
        ST_Transform(boundingbox.geom, 4326),
        polyunion.geom
      ) THEN ST_Difference(
        ST_Transform(boundingbox.geom, 4326),
        polyunion.geom
      )
      ELSE ST_Transform(boundingbox.geom, 4326)
    END AS geom,
    'Out' AS name
  FROM
    boundingbox,
    polyunion

polyunion is the union of the data polygons within a tile represented by its boundingbox. We have two cases: in the first case the boundingbox intersects the data polygons. In this case we take the difference (that is, remove the area of the data from the tile to only have the complement of the data) and assign the ‘Out’ to the field name. In the second case, the two do not intersect, and so we take the full boundingbox and assign ‘Out’ to the field name.

Now we add to this the area where the data is present:

UNION
  SELECT
    ST_Intersection(
      origTable.geom,
      ST_Transform(boundingbox.geom, 4326)
    ) AS geom,
    name
  FROM
    origTable,
    boundingbox

Putting these together and restricting the query to the intersection of the data and the tile, we have:

UNION
  SELECT
    ST_Intersection(
      origTable.geom,
      ST_Transform(boundingbox.geom, 4326)
    ) AS geom,
    name
  FROM
    origTable,
    boundingbox
  SELECT
    CASE
      WHEN ST_Intersects(
        ST_Transform(boundingbox.geom, 4326),
        polyunion.geom
      ) THEN ST_Difference(
        ST_Transform(boundingbox.geom, 4326),
        polyunion.geom
      )
      ELSE ST_Transform(boundingbox.geom, 4326)
    END AS geom,
    'Out' AS name
  FROM
    boundingbox,
    polyunion
  UNION
  SELECT
    ST_Intersection(
      origTable.geom,
      ST_Transform(boundingbox.geom, 4326)
    ) AS geom,
    name
  FROM
    origTable,
    boundingbox
  WHERE
    ST_Intersects(
      origTable.geom,
      ST_Transform(boundingbox.geom, 4326)
    )

That is the unique part. The remainder of the query is the common table expression to create the bounding box using the format required by db2vector:

WITH boundingbox AS(
  SELECT
    ST_MakeEnvelope(
      %(xmin)s,
      %(ymin)s,
      %(xmax)s,
      %(ymax)s,
      3857
    ) AS geom
)

We use two more common table expressions to limit the query to the region of the polygon data and also to collect the polygons as one dissolved polygon to allow taking the difference. One could move these clauses into the main query, but I find it easier to follow the large query when they are separated.

origTable AS (
  SELECT
    ST_Intersection(
      ST_Transform(boundingbox.geom, 4326),
      india_states.geom
    ) AS geom,
    name
  FROM
    demo.india_states,
    boundingbox
  WHERE
    ST_Intersects(
      ST_Transform(boundingbox.geom, 4326),
      india_states.geom
    )
),
polyunion AS (
  SELECT
    ST_Union(origTable.geom) AS geom
  FROM
    origTable,
    boundingbox
  WHERE
    ST_Intersects(
      origTable.geom,
      ST_Transform(boundingbox.geom, 4326)
    )
),

And one last final bit that is specifically for creating vector tiles from PostGIS. If you’re not familiar with vector tiles in this setting, please see the article: https://www.clockworkmicro.com/help/getting-started

mvtgeom AS (
  SELECT
    ST_AsMVTGeom(
      ST_Transform(sourceTable.geom, 3857),
      boundingbox.geom
    ),
    name
  FROM
    sourceTable,
    boundingbox
  WHERE
    ST_Intersects(
      ST_Transform(boundingbox.geom, 4326),
      sourceTable.geom
    )
)
SELECT
  ST_AsMVT(mvtgeom.*)
FROM
  mvtgeom;

Putting all these pieces together, we have:

WITH boundingbox AS(
  SELECT
    ST_MakeEnvelope(
      %(xmin)s,
      %(ymin)s,
      %(xmax)s,
      %(ymax)s,
      3857
    ) AS geom
),
origTable AS (
  SELECT
    ST_Intersection(
      ST_Transform(boundingbox.geom, 4326),
      india_states.geom
    ) AS geom,
    name
  FROM
    demo.india_states,
    boundingbox
  WHERE
    ST_Intersects(
      ST_Transform(boundingbox.geom, 4326),
      india_states.geom
    )
),
polyunion AS (
  SELECT
    ST_Union(origTable.geom) AS geom
  FROM
    origTable,
    boundingbox
  WHERE
    ST_Intersects(
      origTable.geom,
      ST_Transform(boundingbox.geom, 4326)
    )
),
sourceTable AS (
  SELECT
    CASE
      WHEN ST_Intersects(
        ST_Transform(boundingbox.geom, 4326),
        polyunion.geom
      ) THEN ST_Difference(
        ST_Transform(boundingbox.geom, 4326),
        polyunion.geom
      )
      ELSE ST_Transform(boundingbox.geom, 4326)
    END AS geom,
    'Out' AS name
  FROM
    boundingbox,
    polyunion
  UNION
  SELECT
    ST_Intersection(
      origTable.geom,
      ST_Transform(boundingbox.geom, 4326)
    ) AS geom,
    name
  FROM
    origTable,
    boundingbox
  WHERE
    ST_Intersects(
      origTable.geom,
      ST_Transform(boundingbox.geom, 4326)
    )
),
mvtgeom AS (
  SELECT
    ST_AsMVTGeom(
      ST_Transform(sourceTable.geom, 3857),
      boundingbox.geom
    ),
    name
  FROM
    sourceTable,
    boundingbox
  WHERE
    ST_Intersects(
      ST_Transform(boundingbox.geom, 4326),
      sourceTable.geom
    )
)
SELECT
  ST_AsMVT(mvtgeom.*)
FROM
  mvtgeom;

Here is a side-by-side of our final map showing a the states of India and their complement.