So in this post the first in a series of little real world use cases of MapLesotho data, we’re going to be mapping building density in Lesotho.

For those not in the know that means we’re going to split lesotho into little parcels of a fixed size and then count the number of buildings within each parcel.

Assumptions

We will make the following assumptions if you’re following along, that you’ve got:

• a postgres database with postgis loaded.
• an export of Lesotho loaded in with osm2pgsql

The first thing we’re going to do is create a hexgon grid to cover the land area of Lesotho. Although we can do this in QGIS i’m going to show you how to easily create one in PostGIS.

Thankfully CartoDB has already made a super handy function for use to create on of these grid. You should run the code below in your osm database to add the CDB_HexagonGrid() function.

Source Copy from this gist

Now we shall make our hexagon grid table and make one that matches the are of Lesotho.

-- create our hexgrid table
DROP TABLE IF EXISTS hexgrid;
CREATE TABLE hexgrid (
  gid             SERIAL PRIMARY KEY,
  cell           geometry NOT NULL
);

-- add the cells to the grid using lesotho admin boundary
insert into hexgrid(way)
  WITH
 geo as (select way from planet_osm_polygon where osm_id = -2093234),
 grid AS ( SELECT CDB_HexagonGrid(ST_Envelope(geo.way), 2000) AS cell from geo)
SELECT way from grid

And the result:

Building Density per grid cell

So next we’re going to do a little prepartory work to make our queries work a little faster when querying the buildings of Lesotho. There are > 750,000 building in case you didn’t know.

First we will create a new table containing only buildings. You’ll notice we’re transforming our geomtetry, this is to allow us to take advantage of geography queries in later posts. We also get the centre point of the each polygon.

-- get all the buildings into one table
DROP TABLE IF EXISTS buildings;
CREATE TABLE buildings as (
SELECT osm_id, tags, way_area, ST_Transform(way,4326) as way, ST_Centroid(way) as centroid
FROM planet_osm_polygon where building IS NOT NULL
);

Now to help our queries run much faster, we will create some spatial indexs on the geometry and osm_id columns. This may take a few minutes to run through.

-- WE NEED TO MAKE THESE INDEXS SO THINGS RUN FASTER
CREATE INDEX buildings_index
  ON public.buildings
  USING gist
  (way);

CREATE INDEX buildings_centroid_index
  ON public.buildings
  USING gist
  (centroid);

CREATE INDEX buildings_pkey
  ON public.buildings
  USING btree
  (osm_id);

Here’s a quick image of the buildings and their center points.

Now onto the fancy part we’re going to use some spatial joins to count how many centroid points are within each hexagon. Running the query below will add a now column to the hexgrid table for the building count.

ALTER TABLE hexgrid ADD COLUMN agg_building float; --Add a column to the hexagon table
UPDATE hexgrid SET agg_building = 0.0; --Initialize that column to a count of zero
UPDATE hexgrid SET agg_building = temporary_holder.building_count
  FROM
  (
    SELECT
      hexgrid.gid, --Keep the ID for each hexagon to help specify which rows to update
      count( --The aggregation function used in collapsing to just hexagons
          buildings.osm_id --Just a vector of ones meaning yes if building
          ) AS building_count --Give the total a name for easy reference
    FROM
    (
      hexgrid
      INNER JOIN -- Make a new table where each pair of hexagon-gazetteer point is a row
      buildings
      ON ST_Intersects(hexgrid.cell, buildings.centroid) -- Only for pairs that occupy the same space
    )
    GROUP BY hexgrid.gid -- Collapse hexagon-gazetteer point pairs to just hexagons
  ) AS temporary_holder --The placeholder name for the aggregate results
  WHERE hexgrid.gid=temporary_holder.gid --Update rows in the hexagon table with matching id number

Original Query - Rex Douglass

If we bring these layers into QGIS we can then visualise the results:

You can see that the highest densities are surrounding the major settlements and in the north from Maseru to Bera, to Leribe, and Buthe-Buthe.

Have an idea another idea of a #MapLesotho question we could answer with this technique, then let me know on twitter @Rusty1052