Using Uber H3 Indexing Library in Postgres for Geospatial Data Analytics

# Pull any postgis docker image
docker pull postgis/postgis:15-master
# Let's run a demo container where we can install h3 extension
docker run -it --name postgis-h3 -e POSTGRES_PASSWORD=postgres postgis/postgis:15-master

# shell access into the docker
docker exec -it -u root postgis-h3 bash

# Install the dependencies
apt update
apt install -y pip libpq-dev postgresql-server-dev-15
pip install pgxnclient cmake
# Install the extension
pgxn install h3

# Remove the dev dependencies
pip uninstall pgxnclient cmake
apt purge -y libpq-dev postgresql-server-dev-15 pip
exit

docker commit postgis-h3 postgis-h3
docker stop postgis-h3
docker rm postgis-h3
docker run -it --name postgis-h3 -d  -p 6432:5432 -e POSTGRES_PASSWORD=postgres postgis-h3

docker run -it --name postgis-h3 -d -p 6432:5432 -e POSTGRES_PASSWORD=postgres  jsonsingh/postgis-h3

PGPASSWORD=postgres psql -p 6432 -h localhost

-- Create a temporary database
Create database temp;
\c temp;
create extension h3;
create extension h3_postgis CASCADE;
SELECT h3_lat_lng_to_cell(POINT('37.3615593,-122.0553238'), 5);

-- Clean up
\c postgres
Drop database temp

# Just to make our life easy
export PGPASSWORD=postgres
export PGPORT=6432
export PGHOST=localhost

# Needed to load osm data to postgres
apt install osm2pgsql

# our database and extension
psql -c "CREATE database osm_data;"
psql -d osm_data -c "CREATE extension postgis;"

# This will take some time
osm2pgsql -c -d osm_data -x -E 4326 <osm_file>
# Mulitple OSM files
osm2pgsql -a -d osm_data -x -E 4326 <osm_file_2>

PGPASSWORD=postgres psql -p 6432 -h localhost -d osm_data

CREATE TABLE planet_osm_polygon_admin_6 AS SELECT DISTINCT ON (osm_id) * FROM planet_osm_polygon WHERE admin_level = '6';
CREATE TABLE planet_osm_trees AS SELECT * FROM planet_osm_point WHERE "natural" = 'tree';

SELECT COUNT(*) FROM planet_osm_polygon_admin_6;
-- 5518

SELECT COUNT(*) FROM planet_osm_trees;
-- 799130

-- Not let's try to find all the trees in these countries
SELECT
    b.osm_id,
    COUNT(a.osm_id)
FROM
    planet_osm_trees a
    RIGHT JOIN planet_osm_polygon_admin_6 b ON ST_Within (a.way, b.way)
GROUP BY
    b.osm_id
ORDER BY
    COUNT;
-- took 343 seconds

CREATE extension h3;
CREATE extension h3_postgis CASCADE;
-- Let's try the h3 grid. Calculate all the h3 index for all trees
ALTER TABLE planet_osm_trees ADD COLUMN h3_index h3index GENERATED ALWAYS AS (h3_lat_lng_to_cell(way::POINT, 7)) STORED;

-- For polygons we need to do a bit of extra work
-- This function takes a polygon and gives all the h3 index inside that polygon
CREATE OR REPLACE FUNCTION get_h3_indexes(shape geometry, index integer)
  RETURNS h3index[] AS $$
DECLARE
  h3_indexes h3index[];
BEGIN
  SELECT ARRAY(
    SELECT h3_polygon_to_cells(shape, index)
  ) INTO h3_indexes;

  RETURN h3_indexes;
END;
$$ LANGUAGE plpgsql IMMUTABLE;

-- We calculate all the h3 index for all the shapes. This will take a lot of time but it's a one time process
ALTER TABLE planet_osm_polygon_admin_6  ADD COLUMN h3_indexes h3index[] GENERATED ALWAYS AS (get_h3_indexes(way,7)) STORED;
-- 30 seconds

-- Only for visualization purpose
ALTER TABLE planet_osm_trees ADD COLUMN h3_shape geometry GENERATED ALWAYS AS (h3_cell_to_boundary_geometry(h3_lat_lng_to_cell(way::POINT, 7))) STORED;

-- Create a flat table for faster joins and visualize this table in QGIS
CREATE TABLE
    planet_osm_polygon_admin_6_flat AS
SELECT
    h3_index,
    osm_id,
    -- This is just for visualization and not necessary
    h3_cell_to_boundary_geometry (h3_index)
FROM
    (
        SELECT
            UNNEST (h3_indexes) AS h3_index,
            osm_id
        FROM
            planet_osm_polygon_admin_6
    ) AS a;
-- 8 seconds

-- We compare the h3 index and find all the shape ids
SELECT
    b.osm_id,
    COUNT(a.osm_id)
FROM
    planet_osm_trees a
    JOIN planet_osm_polygon_admin_6_flat b ON a.h3_index = b.h3_index
GROUP BY
    b.osm_id
ORDER BY
    COUNT;
-- takes 1 sec

-- create a function which returns compacted h3
CREATE OR REPLACE FUNCTION get_h3_indexes_compact(shape geometry, index integer) RETURNS h3index[] AS $$ DECLARE h3_indexes h3index[];
BEGIN WITH cells AS (
  SELECT
    h3_polygon_to_cells(shape, index) AS cell_array
)
SELECT
  Array(
    SELECT
      h3_compact_cells(
        array_agg(cell_array)
      )
    FROM
      cells
  ) into h3_indexes;
RETURN h3_indexes;
END;
$$ LANGUAGE plpgsql IMMUTABLE;


-- Let's create the compacted array indice
ALTER TABLE
  planet_osm_polygon_admin_6
ADD
  COLUMN h3_indexes_compacted h3index[] GENERATED ALWAYS AS (
    get_h3_indexes_compact(way, 7)
  ) STORED;


-- Visualize the below table in QGIS
CREATE TABLE
    planet_osm_polygon_admin_6_flat_compact AS
SELECT
    h3_index,
    osm_id,
    h3_cell_to_boundary_geometry (h3_index)
FROM
    (
        SELECT
            UNNEST (h3_indexes_compacted) AS h3_index,
            osm_id
        FROM
            planet_osm_polygon_admin_6
    ) AS a;