diff --git a/app/airq/models/relations.py b/app/airq/models/relations.py index 3422281..456831b 100644 --- a/app/airq/models/relations.py +++ b/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: diff --git a/app/airq/models/sensors.py b/app/airq/models/sensors.py index 644ef35..a7e0a86 100644 --- a/app/airq/models/sensors.py +++ b/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"" diff --git a/app/airq/sync/purpleair.py b/app/airq/sync/purpleair/__init__.py similarity index 77% rename from app/airq/sync/purpleair.py rename to app/airq/sync/purpleair/__init__.py index 73437aa..164bfe5 100644 --- a/app/airq/sync/purpleair.py +++ b/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() diff --git a/app/airq/sync/purpleair/metrics.py b/app/airq/sync/purpleair/metrics.py new file mode 100644 index 0000000..9a0e164 --- /dev/null +++ b/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() diff --git a/app/migrations/versions/28fc973c085a_add_coords_to_sensors.py b/app/migrations/versions/28fc973c085a_add_coords_to_sensors.py new file mode 100644 index 0000000..784d9f3 --- /dev/null +++ b/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 ### diff --git a/app/tests/test_clients.py b/app/tests/test_clients.py index 7693bed..549c11e 100644 --- a/app/tests/test_clients.py +++ b/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 ) diff --git a/app/tests/test_sms.py b/app/tests/test_sms.py index 83b1d0e..60379df 100644 --- a/app/tests/test_sms.py +++ b/app/tests/test_sms.py @@ -29,7 +29,7 @@ def test_get_quality(self): self.assertEqual(1, Client.query.count()) self.assertEqual(1, Event.query.count()) self.assert_twilio_response( - "Welcome to Hazebot! We'll send you alerts when air quality in Portland 97204 changes category. Air quality is now GOOD (AQI 42).\n" + "Welcome to Hazebot! We'll send you alerts when air quality in Portland 97204 changes category. Air quality is now GOOD (AQI 44).\n" "\n" 'Save this contact and text us your zipcode whenever you\'d like an instant update. And you can always text "M" to see the whole menu.', response.data, @@ -37,7 +37,7 @@ def test_get_quality(self): ) client_id = Client.query.filter_by(identifier="+12222222222").first().id - self.assert_event(client_id, EventType.QUALITY, zipcode="97204", pm25=9.875) + self.assert_event(client_id, EventType.QUALITY, zipcode="97204", pm25=10.557) response = self.client.post( "/sms/en", data={"Body": "2", "From": "+12222222222"} @@ -46,12 +46,12 @@ def test_get_quality(self): self.assertEqual(1, Client.query.count()) self.assertEqual(2, Event.query.count()) self.assert_twilio_response( - "Portland 97204 is GOOD (AQI 42).\n" + "Portland 97204 is GOOD (AQI 44).\n" "\n" 'Text "M" for Menu, "E" to end alerts.', response.data, ) - self.assert_event(client_id, EventType.LAST, zipcode="97204", pm25=9.875) + self.assert_event(client_id, EventType.LAST, zipcode="97204", pm25=10.557) response = self.client.post( "/sms/en", data={"Body": "1", "From": "+12222222222"} @@ -62,14 +62,14 @@ def test_get_quality(self): self.assert_twilio_response( "GOOD (AQI: 0 - 50) means air quality is considered satisfactory, and air pollution poses little or no risk.\n" "\n" - "Average PM2.5 from 8 sensors near 97204 is 9.875 ug/m^3.", + "Average PM2.5 from 8 sensors near 97204 is 10.557 ug/m^3.", response.data, ) self.assert_event( client_id, EventType.DETAILS, zipcode="97204", - pm25=9.875, + pm25=10.557, num_sensors=8, recommendations=[], ) @@ -81,7 +81,7 @@ def test_get_quality(self): self.assertEqual(2, Client.query.count()) self.assertEqual(4, Event.query.count()) self.assert_twilio_response( - "Welcome to Hazebot! We'll send you alerts when air quality in Molalla 97038 changes category. Air quality is now MODERATE (AQI 74).\n" + "Welcome to Hazebot! We'll send you alerts when air quality in Molalla 97038 changes category. Air quality is now MODERATE (AQI 69).\n" "\n" 'Save this contact and text us your zipcode whenever you\'d like an instant update. And you can always text "M" to see the whole menu.', response.data, @@ -89,7 +89,7 @@ def test_get_quality(self): ) client_id = Client.query.filter_by(identifier="+13333333333").first().id - self.assert_event(client_id, EventType.QUALITY, zipcode="97038", pm25=22.9) + self.assert_event(client_id, EventType.QUALITY, zipcode="97038", pm25=20.416) response = self.client.post( "/sms/en", data={"Body": "2", "From": "+13333333333"} @@ -98,12 +98,12 @@ def test_get_quality(self): self.assertEqual(2, Client.query.count()) self.assertEqual(5, Event.query.count()) self.assert_twilio_response( - "Molalla 97038 is MODERATE (AQI 74).\n" + "Molalla 97038 is MODERATE (AQI 69).\n" "\n" 'Text "M" for Menu, "E" to end alerts.', response.data, ) - self.assert_event(client_id, EventType.LAST, zipcode="97038", pm25=22.9) + self.assert_event(client_id, EventType.LAST, zipcode="97038", pm25=20.416) response = self.client.post( "/sms/en", data={"Body": "1", "From": "+13333333333"} @@ -117,19 +117,19 @@ def test_get_quality(self): "Here are the closest places with better air quality:" "\n" " - Estacada 97023: GOOD (16.7 mi)\n" + " - West Linn 97068: GOOD (17.3 mi)\n" " - Gladstone 97027: GOOD (18.5 mi)\n" - " - Eagle Creek 97022: GOOD (20.0 mi)\n" "\n" - "Average PM2.5 from 3 sensors near 97038 is 22.9 ug/m^3.", + "Average PM2.5 from 3 sensors near 97038 is 20.416 ug/m^3.", response.data, ) self.assert_event( client_id, EventType.DETAILS, zipcode="97038", - pm25=22.9, + pm25=20.416, num_sensors=3, - recommendations=["97023", "97027", "97022"], + recommendations=["97023", "97068", "97027"], ) self.clock.advance() @@ -140,12 +140,12 @@ def test_get_quality(self): self.assertEqual(2, Client.query.count()) self.assertEqual(7, Event.query.count()) self.assert_twilio_response( - "Molalla 97038 is MODERATE (AQI 74).\n" + "Molalla 97038 is MODERATE (AQI 69).\n" "\n" 'Text "M" for Menu, "E" to end alerts.', response.data, ) - self.assert_event(client_id, EventType.QUALITY, zipcode="97038", pm25=22.9) + self.assert_event(client_id, EventType.QUALITY, zipcode="97038", pm25=20.416) self.clock.advance() response = self.client.post( @@ -155,12 +155,12 @@ def test_get_quality(self): self.assertEqual(2, Client.query.count()) self.assertEqual(8, Event.query.count()) self.assert_twilio_response( - "Portland 97204 is GOOD (AQI 42).\n" + "Portland 97204 is GOOD (AQI 44).\n" "\n" "You are now receiving alerts for 97204.", response.data, ) - self.assert_event(client_id, EventType.QUALITY, zipcode="97204", pm25=9.875) + self.assert_event(client_id, EventType.QUALITY, zipcode="97204", pm25=10.557) def test_get_menu(self): expected_response = ( @@ -259,7 +259,7 @@ def test_unsubscribe(self): self.assertEqual(1, Client.query.count()) self.assertEqual(1, Event.query.count()) self.assert_twilio_response( - "Welcome to Hazebot! We'll send you alerts when air quality in Portland 97204 changes category. Air quality is now GOOD (AQI 42).\n" + "Welcome to Hazebot! We'll send you alerts when air quality in Portland 97204 changes category. Air quality is now GOOD (AQI 44).\n" "\n" 'Save this contact and text us your zipcode whenever you\'d like an instant update. And you can always text "M" to see the whole menu.', response.data, @@ -269,7 +269,7 @@ def test_unsubscribe(self): client = Client.query.first() self.assertEqual("97204", client.zipcode.zipcode) self.assertEqual(0, client.alerts_disabled_at) - self.assert_event(client.id, EventType.QUALITY, zipcode="97204", pm25=9.875) + self.assert_event(client.id, EventType.QUALITY, zipcode="97204", pm25=10.557) alerts_disabled_at = self.clock.advance().timestamp() response = self.client.post( @@ -321,7 +321,7 @@ def test_unsubscribe(self): self.assertEqual(1, Client.query.count()) self.assertEqual(4, Event.query.count()) self.assert_twilio_response( - "Portland 97204 is GOOD (AQI 42).\n" + "Portland 97204 is GOOD (AQI 44).\n" "\n" 'Alerting is disabled. Text "Y" to re-enable alerts when air quality changes.', response.data, @@ -330,7 +330,7 @@ def test_unsubscribe(self): client = Client.query.first() self.assertEqual("97204", client.zipcode.zipcode) self.assertEqual(alerts_disabled_at, client.alerts_disabled_at) - self.assert_event(client.id, EventType.QUALITY, zipcode="97204", pm25=9.875) + self.assert_event(client.id, EventType.QUALITY, zipcode="97204", pm25=10.557) self.clock.advance() response = self.client.post( @@ -565,7 +565,7 @@ def test_translations(self): self.assertEqual(1, Client.query.count()) self.assertEqual(1, Event.query.count()) self.assert_twilio_response( - "¡Bienvenido a Hazebot! Le enviaremos avisos cuando la calidad del aire en Portland 97204 cambie de categoría. La calidad del aire ahora es BUENO (AQI 42).\n" + "¡Bienvenido a Hazebot! Le enviaremos avisos cuando la calidad del aire en Portland 97204 cambie de categoría. La calidad del aire ahora es BUENO (AQI 44).\n" "\n" 'Guardar este contacto y enviarnos un mensaje de texto con su código postal cuando desee una actualización instantánea. Y siempre puede enviar un mensaje de texto con "M" para ver el menú completo.', response.data, diff --git a/app/tests/test_zipcodes.py b/app/tests/test_zipcodes.py index 5fa2375..fd883e8 100644 --- a/app/tests/test_zipcodes.py +++ b/app/tests/test_zipcodes.py @@ -45,8 +45,8 @@ def test_get_recommendations(self): self.assertListEqual( [ Zipcode.query.filter_by(zipcode="97023").first(), + Zipcode.query.filter_by(zipcode="97068").first(), Zipcode.query.filter_by(zipcode="97027").first(), - Zipcode.query.filter_by(zipcode="97022").first(), ], zipcode.get_recommendations(3, ConversionFactor.NONE), )