Draft: Use Voronoi Diagram for Better AQI Estimates

app/airq/models/relations.py

@@ -10,6 +10,7 @@ class SensorZipcodeRelation(db.Model):  # type: ignore
     zipcode_id = db.Column(
         db.Integer(), db.ForeignKey("zipcodes.id"), nullable=False, primary_key=True
     )
+    # Unit: KM
     distance = db.Column(db.Float(), nullable=False)
 
     def __repr__(self) -> str:

app/airq/models/sensors.py

@@ -1,4 +1,5 @@
 from flask_sqlalchemy import BaseQuery
+from geoalchemy2 import Geometry
 
 from airq.config import db
 
@@ -27,6 +28,8 @@ class Sensor(db.Model):  # type: ignore
     updated_at = db.Column(db.Integer(), nullable=False)
     latitude = db.Column(db.Float(), nullable=False)
     longitude = db.Column(db.Float(), nullable=False)
+
+    # TODO: We can probably use the "coordinates" field instead of geohashing.
     geohash_bit_1 = db.Column(db.String(), nullable=False)
     geohash_bit_2 = db.Column(db.String(), nullable=False)
     geohash_bit_3 = db.Column(db.String(), nullable=False)
@@ -40,6 +43,8 @@ class Sensor(db.Model):  # type: ignore
     geohash_bit_11 = db.Column(db.String(), nullable=False)
     geohash_bit_12 = db.Column(db.String(), nullable=False)
 
+    coordinates = db.Column(Geometry("POINT"), nullable=True)
+
     def __repr__(self) -> str:
         return f"<Sensor {self.id}: {self.latest_reading}>"
 

app/airq/sync/purpleair.py → app/airq/sync/purpleair/__init__.py

@@ -1,4 +1,3 @@
-import collections
 import geohash
 import json
 import logging
@@ -22,6 +21,7 @@
 from airq.models.relations import SensorZipcodeRelation
 from airq.models.sensors import Sensor
 from airq.models.zipcodes import Zipcode
+from airq.sync.purpleair import metrics
 
 
 # Try to get at least 8 readings per zipcode.
@@ -166,9 +166,13 @@ def _sensors_sync(
                 data.update(
                     latitude=latitude,
                     longitude=longitude,
+                    coordinates=f"POINT({longitude} {latitude})",
                     **{f"geohash_bit_{i}": c for i, c in enumerate(gh, start=1)},
                 )
                 moved_sensor_ids.append(result["sensor_index"])
+            elif sensor.coordinates is None:
+                # Sensor wasn't moved, but we haven't filled in its coordinates field yet.
+                data["coordinates"] = f"POINT({longitude} {latitude})"
 
             if sensor:
                 updates.append(data)
@@ -250,82 +254,6 @@ def _relations_sync(moved_sensor_ids: typing.List[int]):
             db.session.commit()
 
 
-def _metrics_sync():
-    logger = get_celery_logger()
-    updates = []
-    ts = now()
-
-    zipcodes_to_sensors = collections.defaultdict(list)
-    for zipcode_id, latest_reading, humidity, pm_cf_1, sensor_id, distance in (
-        Sensor.query.join(SensorZipcodeRelation)
-        .filter(Sensor.updated_at > ts.timestamp() - (30 * 60))
-        .with_entities(
-            SensorZipcodeRelation.zipcode_id,
-            Sensor.latest_reading,
-            Sensor.humidity,
-            Sensor.pm_cf_1,
-            Sensor.id,
-            SensorZipcodeRelation.distance,
-        )
-        .all()
-    ):
-        zipcodes_to_sensors[zipcode_id].append(
-            (latest_reading, humidity, pm_cf_1, sensor_id, distance)
-        )
-
-    for zipcode_id, sensor_tuples in zipcodes_to_sensors.items():
-        pm_25_readings: typing.List[float] = []
-        pm_cf_1_readings: typing.List[float] = []
-        humidities: typing.List[float] = []
-        closest_reading = float("inf")
-        farthest_reading = 0.0
-        sensor_ids: typing.List[int] = []
-        for pm_25, humidity, pm_cf_1, sensor_id, distance in sorted(
-            sensor_tuples, key=lambda s: s[-1]
-        ):
-            if (
-                len(pm_25_readings) < DESIRED_NUM_READINGS
-                or distance < DESIRED_READING_DISTANCE_KM
-            ):
-                pm_25_readings.append(pm_25)
-                humidities.append(humidity)
-                pm_cf_1_readings.append(pm_cf_1)
-                sensor_ids.append(sensor_id)
-                closest_reading = min(distance, closest_reading)
-                farthest_reading = max(distance, farthest_reading)
-            else:
-                break
-
-        if pm_25_readings:
-            num_sensors = len(pm_25_readings)
-            pm25 = round(sum(pm_25_readings) / num_sensors, ndigits=3)
-            humidity = round(sum(humidities) / num_sensors, ndigits=3)
-            pm_cf_1 = round(sum(pm_cf_1_readings) / num_sensors, ndigits=3)
-            min_sensor_distance = round(closest_reading, ndigits=3)
-            max_sensor_distance = round(farthest_reading, ndigits=3)
-            details = {
-                "num_sensors": num_sensors,
-                "min_sensor_distance": min_sensor_distance,
-                "max_sensor_distance": max_sensor_distance,
-                "sensor_ids": sensor_ids,
-            }
-            updates.append(
-                {
-                    "id": zipcode_id,
-                    "pm25": pm25,
-                    "humidity": humidity,
-                    "pm_cf_1": pm_cf_1,
-                    "pm25_updated_at": ts.timestamp(),
-                    "metrics_data": details,
-                }
-            )
-
-    logger.info("Updating %s zipcodes", len(updates))
-    for mappings in chunk_list(updates, batch_size=5000):
-        db.session.bulk_update_mappings(Zipcode, mappings)
-        db.session.commit()
-
-
 def _send_alerts():
     logger = get_celery_logger()
     num_sent = 0
@@ -372,7 +300,7 @@ def purpleair_sync():
         _relations_sync(moved_sensor_ids)
 
     logger.info("Syncing metrics")
-    _metrics_sync()
+    metrics.update()
 
     logger.info("Sending alerts")
     _send_alerts()

app/airq/sync/purpleair/metrics.py

@@ -0,0 +1,180 @@
+import textwrap
+import time
+
+from sqlalchemy import text
+
+from airq.celery import get_celery_logger
+from airq.config import db
+from airq.lib.clock import now
+from airq.lib.util import chunk_list
+from airq.models.zipcodes import Zipcode
+
+# Try to get at least 8 readings per zipcode.
+DESIRED_NUM_READINGS = 8
+
+# Allow any number of readings within 2.5km from the zipcode centroid.
+DESIRED_READING_DISTANCE_KM = 2.5
+
+
+def _execute_query(timestamp):
+    # This approach is really slow (about 5 mins) and hard to test
+    # since it's so complex. How can we make this better?
+    #
+    # For testability: it would be useful to test this in a single
+    # zipcode containing a few sensors. That would allow us to assert
+    # that the algorithm works correctly:
+    # - When there are two equidistant sensors, it assigns equal weight to each.
+    # - When there are two sensors which are not equidistant, the one closer
+    #   to the origin gets more weight.
+    # - When there are three equidistant sensors, it assigns equal weight to each.
+    # - When there are three sensors which are not equidistant, the closer
+    #   sensors are assigned more weight.
+    #
+    # For speed: ¯\_(ツ)_/¯
+
+    sql = text(
+        textwrap.dedent(
+            """
+    -- Compute the Voronoi Diagram for the set of sensors.
+    WITH voronoi_cells AS (
+        SELECT
+            ST_Intersection(
+                (ST_Dump(
+                    ST_VoronoiPolygons(
+                        ST_Collect(coordinates)
+                    ) 
+                )).geom,
+                ST_MakeEnvelope(-180, -90, 180, 90)
+            ) as cell
+        FROM sensors
+        WHERE coordinates IS NOT NULL
+        AND updated_at > :updated_at
+    ),
+
+    -- Map each Voronoi cell to the sensor it contains.
+    -- This is actually the slowest part of this query
+    -- but I'm not sure how to speed it up.
+    sensors_with_cells AS (
+        SELECT 
+            s.id,
+            s.latest_reading,
+            s.humidity,
+            s.pm_cf_1,
+            ST_Area(v.cell) as area
+        FROM sensors s
+        JOIN voronoi_cells v
+        ON ST_Within(s.coordinates, v.cell)
+    ),
+
+    -- Find the distance of the eighth closest sensor to each zipcode.
+    -- We use this to ensure that if a zipcode has no sensors within
+    -- a 2.5 KM radius, we can search outside that radius for at most
+    -- eight sensors. This won't result in us choosing sensors really
+    -- far away because the `sensors_zipcodes` only contains relations
+    -- between sensors and zipcodes at most 20 KM apart.
+    zipcodes_to_distance AS (
+        SELECT
+            z.id,
+            (
+                SELECT MAX(s2.distance)
+                FROM (
+                    SELECT s1.distance
+                    FROM sensors_zipcodes s1
+                    WHERE s1.zipcode_id = z.id
+                    ORDER BY s1.distance
+                    LIMIT :desired_num_readings
+                ) s2
+            ) as distance_to_eighth_closest_sensor
+        FROM zipcodes z
+        GROUP BY z.id
+    ) 
+
+    -- For each zipcode, compute metrics for all eligible sensors.
+    SELECT
+        zd.id,
+        SUM(sc.latest_reading * sc.area) / SUM(sc.area) as pm25,
+        SUM(sc.humidity * sc.area) / SUM(sc.area) as humidity,
+        SUM(sc.pm_cf_1 * sc.area) / SUM(sc.area) as pm_cf_1,
+        MAX(sz.distance) AS max_distance,
+        MIN(sz.distance) AS min_distance,
+        COUNT(sc.id) AS num_sensors,
+        ARRAY_AGG(sc.id) AS sensor_ids
+    FROM zipcodes_to_distance zd
+    JOIN sensors_zipcodes sz
+    ON sz.zipcode_id = zd.id
+    JOIN sensors_with_cells sc
+    ON sc.id = sz.sensor_id
+
+    -- Include all sensors within 2.5 KM of the zipcode's
+    -- centroid. If there are fewer than 8 sensors within
+    -- 2.5 KM, include the closest eight sensors.
+    WHERE sz.distance <= GREATEST(zd.distance_to_eighth_closest_sensor, :desired_reading_distance_km)
+
+    GROUP BY zd.id
+    """
+        )
+    )
+
+    return db.engine.execute(
+        sql,
+        {
+            "desired_num_readings": DESIRED_NUM_READINGS,
+            "desired_reading_distance_km": DESIRED_READING_DISTANCE_KM,
+            "updated_at": timestamp - (30 * 60),
+        },
+    )
+
+
+def _compute_updates():
+    logger = get_celery_logger()
+
+    ts = now()
+    start_ts = time.perf_counter()
+    rows = _execute_query(ts.timestamp())
+    end_ts = time.perf_counter()
+    duration = end_ts - start_ts
+    logger.info("executed sql in %f seconds", duration)
+
+    updates = []
+    for row in rows:
+        (
+            zipcode_id,
+            pm25,
+            humidity,
+            pm_cf_1,
+            max_sensor_distance,
+            min_sensor_distance,
+            num_sensors,
+            sensor_ids,
+        ) = row
+
+        details = {
+            "num_sensors": num_sensors,
+            "min_sensor_distance": min_sensor_distance,
+            "max_sensor_distance": max_sensor_distance,
+            "sensor_ids": sensor_ids,
+        }
+
+        updates.append(
+            {
+                "id": zipcode_id,
+                "pm25": round(pm25, ndigits=3),
+                "humidity": round(humidity, ndigits=3),
+                "pm_cf_1": round(pm_cf_1, ndigits=3),
+                "pm25_updated_at": ts.timestamp(),
+                "metrics_data": details,
+            }
+        )
+
+    return updates
+
+
+def update():
+    updates = _compute_updates()
+
+    logger = get_celery_logger()
+    logger.info("Updating %d zipcodes", len(updates))
+
+    for mappings in chunk_list(updates, batch_size=5000):
+        db.session.bulk_update_mappings(Zipcode, mappings)
+        db.session.commit()

app/migrations/versions/28fc973c085a_add_coords_to_sensors.py

@@ -0,0 +1,38 @@
+"""add coords to sensors
+
+Revision ID: 28fc973c085a
+Revises: 5f3e5ff4f100
+Create Date: 2022-02-13 07:04:14.433225
+
+"""
+from alembic import op
+import sqlalchemy as sa
+from geoalchemy2.types import Geometry
+
+
+# revision identifiers, used by Alembic.
+revision = "28fc973c085a"
+down_revision = "5f3e5ff4f100"
+branch_labels = None
+depends_on = None
+
+
+def upgrade():
+    # ### commands auto generated by Alembic - please adjust! ###
+    op.add_column(
+        "sensors",
+        sa.Column(
+            "coordinates",
+            Geometry(
+                geometry_type="POINT", from_text="ST_GeomFromEWKT", name="geometry"
+            ),
+            nullable=True,
+        ),
+    )
+    # ### end Alembic commands ###
+
+
+def downgrade():
+    # ### commands auto generated by Alembic - please adjust! ###
+    op.drop_column("sensors", "coordinates")
+    # ### end Alembic commands ###

app/tests/test_clients.py

@@ -212,7 +212,7 @@ def test_disable_alerts(self):
         client = self._make_client(last_pm25=self.zipcode.pm25)
         client.disable_alerts()
         self.assertEqual(self.timestamp, client.alerts_disabled_at)
-        self.assertEqual(9.875, client.last_pm25)
+        self.assertEqual(10.557, client.last_pm25)
         self.assert_event(
             client.id, EventType.UNSUBSCRIBE, zipcode=client.zipcode.zipcode
         )