From 36c6bd7e4662be4a7ce62a1f4aad44cab95c2cca Mon Sep 17 00:00:00 2001 From: LongHairedHacker Date: Sun, 16 Jan 2022 16:10:37 +0100 Subject: [PATCH] Initial commit --- pressure_wave_data.ipynb | 409 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 409 insertions(+) create mode 100644 pressure_wave_data.ipynb diff --git a/pressure_wave_data.ipynb b/pressure_wave_data.ipynb new file mode 100644 index 0000000..de64f4b --- /dev/null +++ b/pressure_wave_data.ipynb @@ -0,0 +1,409 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "4dbd5465", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Defaulting to user installation because normal site-packages is not writeable\n", + "Requirement already satisfied: geopy in /home/sebastian/.local/lib/python3.10/site-packages (2.2.0)\n", + "Requirement already satisfied: geographiclib<2,>=1.49 in /home/sebastian/.local/lib/python3.10/site-packages (from geopy) (1.52)\n", + "Defaulting to user installation because normal site-packages is not writeable\n", + "Requirement already satisfied: ipyleaflet in /home/sebastian/.local/lib/python3.10/site-packages (0.15.0)\n", + "Requirement already satisfied: traittypes<3,>=0.2.1 in /home/sebastian/.local/lib/python3.10/site-packages (from ipyleaflet) (0.2.1)\n", + "Requirement already satisfied: ipywidgets<8,>=7.6.0 in /usr/lib/python3.10/site-packages (from ipyleaflet) (7.6.5)\n", + "Requirement already satisfied: xyzservices>=2021.8.1 in /home/sebastian/.local/lib/python3.10/site-packages (from ipyleaflet) (2021.11.0)\n", + "Requirement already satisfied: traitlets>=4.2.2 in /usr/lib/python3.10/site-packages (from traittypes<3,>=0.2.1->ipyleaflet) (5.1.0)\n", + "Defaulting to user installation because normal site-packages is not writeable\n", + "Collecting tabulate\n", + " Downloading tabulate-0.8.9-py3-none-any.whl (25 kB)\n", + "Installing collected packages: tabulate\n", + "\u001b[33m WARNING: The script tabulate is installed in '/home/sebastian/.local/bin' which is not on PATH.\n", + " Consider adding this directory to PATH or, if you prefer to suppress this warning, use --no-warn-script-location.\u001b[0m\n", + "Successfully installed tabulate-0.8.9\n" + ] + } + ], + "source": [ + "import sys\n", + "!{sys.executable} -m pip install geopy\n", + "!{sys.executable} -m pip install ipyleaflet\n", + "!{sys.executable} -m pip install tabulate" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "1ced7adc", + "metadata": {}, + "outputs": [], + "source": [ + "from datetime import datetime\n", + "from pytz import timezone\n", + "\n", + "data_points = [\n", + " (\"Kaiserslautern\", datetime(2022, 1, 16, 20, 30, 0, tzinfo=timezone(\"Europe/Berlin\")), \"https://chaos.social/@sebastian/107628740679869429\"),\n", + " (\"Aachen\", datetime(2022, 1, 16, 20, 25, 0, tzinfo=timezone(\"Europe/Berlin\")), \"https://chaos.social/@trilader/107628784664752716\"),\n", + " (\"Hamburg\", datetime(2022, 1, 16, 20, 10, 0, tzinfo=timezone(\"Europe/Berlin\")),\"https://chaos.social/@tsia_/107628923762308800\"),\n", + " (\"Edinburgh\", datetime(2022, 1, 16, 18, 45, 0, tzinfo=timezone(\"UTC\")), \"https://toot.cat/@river/107629146900457261\"),\n", + " (\"Copenhagen\", datetime(2022, 1, 16, 19, 55, 0, tzinfo=timezone(\"CET\")), \"https://mcd.dk/@rune/107629163090789143\"),\n", + " (\"Augsburg\", datetime(2022, 1, 16, 20, 38, 0, tzinfo=timezone(\"Europe/Berlin\")), \"https://chaos.social/@phjl/107629214329242643\"),\n", + " (\"Frankfurt\", datetime(2022, 1, 16, 20, 25, 0, tzinfo=timezone(\"Europe/Berlin\")), \"https://mastodon.social/@hko/107629722849924036\"),\n", + " (\"Melbourne\", datetime(2022, 1, 16, 8, 0, 0, tzinfo=timezone(\"UTC\")), \"https://aus.social/@futzle/107631253370573653\"),\n", + "]" + ] + }, + { + "cell_type": "markdown", + "id": "1f554c69", + "metadata": {}, + "source": [ + "Since we know roughly where and when the explosion happened we can do a quick sanity check of our data points." + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "id": "6b525fc1", + "metadata": {}, + "outputs": [], + "source": [ + "from geopy.distance import distance\n", + "from geopy.geocoders import Nominatim\n", + "\n", + "geolocator = Nominatim(user_agent=\"pressure_wave_triangulation\")\n", + "\n", + "ORIGIN = (-20.5, -175.4)\n", + "T_ZERO = datetime(2022, 1, 16, 4, 27, 0, tzinfo=timezone(\"UTC\"))\n", + "\n", + "SPEED_OF_SOUND = 343.0 # Speed of sound in air at 20°C\n", + "\n", + "enhanced_data_points = []\n", + "\n", + "for location, ts, url in data_points:\n", + " geolocation = geolocator.geocode(location)\n", + " coordinates = (geolocation.latitude, geolocation.longitude)\n", + " dist = distance(ORIGIN, coordinates)\n", + " \n", + " delta_t = ts - T_ZERO\n", + " \n", + " speed = dist.meters / delta_t.total_seconds()\n", + " mach = speed / SPEED_OF_SOUND\n", + "\n", + " enhanced_data_points += [(location, coordinates, ts, dist, delta_t, speed, mach, url)]\n" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "id": "842f682a", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
Location Coordinates Timestamp Distance Time differenceSpeed [m/s] Mach Number Toot
Kaiserslautern49.443217 7.768995 2022-01-16 20:30:00+00:5316780.572549711476 km15:10:00 307.3364935844593 0.8960247626369076https://chaos.social/@sebastian/107628740679869429
Aachen 50.776351 6.083862 2022-01-16 20:25:00+00:5316641.951165318445 km15:05:00 306.4816052544834 0.8935323768352286https://chaos.social/@trilader/107628784664752716
Hamburg 53.550341 10.000654 2022-01-16 20:10:00+00:5316307.144415681687 km14:50:00 305.3772362487207 0.8903126421245501https://chaos.social/@tsia_/107628923762308800
Edinburgh 55.953346 -3.188375 2022-01-16 18:45:00+00:0016015.614674425336 km14:18:00 311.1036261543383 0.9070076564266423https://toot.cat/@river/107629146900457261
Copenhagen 55.686724 12.570072 2022-01-16 19:55:00+01:0016042.101561186751 km14:28:00 308.028063770866940.8980410022474254https://mcd.dk/@rune/107629163090789143
Augsburg 48.366804 10.898697 2022-01-16 20:38:00+00:5316861.97644995329 km 15:18:00 306.136101124787440.8925250761655611https://chaos.social/@phjl/107629214329242643
Frankfurt 50.110644 8.682092 2022-01-16 20:25:00+00:5316699.02015177992 km 15:05:00 307.532599480293240.8965964999425459https://mastodon.social/@hko/107629722849924036
Melbourne -37.814218 144.9631612022-01-16 08:00:00+00:004264.645211600706 km 3:33:00 333.6968084194606 0.9728769924765615https://aus.social/@futzle/107631253370573653
" + ], + "text/plain": [ + "'\\n\\n\\n\\n\\n\\n\\n\\n\\n\\n\\n\\n
Location Coordinates Timestamp Distance Time differenceSpeed [m/s] Mach Number Toot
Kaiserslautern49.443217 7.768995 2022-01-16 20:30:00+00:5316780.572549711476 km15:10:00 307.3364935844593 0.8960247626369076https://chaos.social/@sebastian/107628740679869429
Aachen 50.776351 6.083862 2022-01-16 20:25:00+00:5316641.951165318445 km15:05:00 306.4816052544834 0.8935323768352286https://chaos.social/@trilader/107628784664752716
Hamburg 53.550341 10.000654 2022-01-16 20:10:00+00:5316307.144415681687 km14:50:00 305.3772362487207 0.8903126421245501https://chaos.social/@tsia_/107628923762308800
Edinburgh 55.953346 -3.188375 2022-01-16 18:45:00+00:0016015.614674425336 km14:18:00 311.1036261543383 0.9070076564266423https://toot.cat/@river/107629146900457261
Copenhagen 55.686724 12.570072 2022-01-16 19:55:00+01:0016042.101561186751 km14:28:00 308.028063770866940.8980410022474254https://mcd.dk/@rune/107629163090789143
Augsburg 48.366804 10.898697 2022-01-16 20:38:00+00:5316861.97644995329 km 15:18:00 306.136101124787440.8925250761655611https://chaos.social/@phjl/107629214329242643
Frankfurt 50.110644 8.682092 2022-01-16 20:25:00+00:5316699.02015177992 km 15:05:00 307.532599480293240.8965964999425459https://mastodon.social/@hko/107629722849924036
Melbourne -37.814218 144.9631612022-01-16 08:00:00+00:004264.645211600706 km 3:33:00 333.6968084194606 0.9728769924765615https://aus.social/@futzle/107631253370573653
'" + ] + }, + "execution_count": 34, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import tabulate\n", + "\n", + "table_data = [(\"Location\", \"Coordinates\", \"Timestamp\", \"Distance\", \"Time difference\", \"Speed [m/s]\", \"Mach Number\", \"Toot\")]\n", + "for location, coordinates, ts, dist, delta_t, speed, mach, url in enhanced_data_points:\n", + " table_data += [(location, \"%f %f\" % coordinates, ts, dist, delta_t, speed, mach, url)]\n", + "\n", + "table = tabulate.tabulate(table_data, tablefmt='html')\n", + "table\n" + ] + }, + { + "cell_type": "markdown", + "id": "5341fa32", + "metadata": {}, + "source": [ + "That looks quite consistent, considering the we pretty much guestimated the timestamps from screenshots." + ] + }, + { + "cell_type": "code", + "execution_count": 120, + "id": "5ff0ee40", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Kaiserslautern 322.0800766318187\n", + "Aachen 306.3924371618201\n", + "Hamburg 282.48002118519315\n", + "Edinburgh 315.16138946596874\n", + "Copenhagen 283.03327184485505\n", + "Augsburg 314.29962074367353\n", + "Frankfurt 319.0927106279109\n" + ] + } + ], + "source": [ + "relative_speed_data = []\n", + "\n", + "for location, coordinates, ts, _, _, _, _, url in enhanced_data_points:\n", + " count = 0\n", + " speed = 0.0\n", + " for _, other_coordinates, other_ts, _, _, _, _, _ in enhanced_data_points:\n", + " relative_delay = other_ts - ts\n", + " if abs(relative_delay.total_seconds()) < 1:\n", + " continue\n", + " \n", + " relative_dist = distance(coordinates, other_coordinates)\n", + " relative_speed = relative_dist.meters / abs(relative_delay.total_seconds())\n", + " \n", + " if relative_speed > 350:\n", + " continue\n", + " \n", + " speed += relative_speed\n", + " count += 1.0\n", + " \n", + " \n", + " if count < 1.0:\n", + " continue \n", + " \n", + " speed = speed / count\n", + " print(location, speed)\n", + " \n", + " relative_speed_data += [(location, coordinates, ts, speed)]\n", + " \n", + " " + ] + }, + { + "cell_type": "markdown", + "id": "5d40d782", + "metadata": {}, + "source": [ + "Some cleanup was required, so I ended up dropping everthing that is much faster then the speed of sound in air.\n", + "That should give use something to work with.\n", + "\n", + "Also the leaftlet map does something weird with the scaling.\n", + "If we cicle each of our locations with the distance we got from geopy, the circles never really intersect." + ] + }, + { + "cell_type": "code", + "execution_count": 121, + "id": "1a43cbc5", + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "175f37ff4cf048d0b886861621cb4598", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Map(center=[-37.8142176, 144.9631608], controls=(ZoomControl(options=['position', 'zoom_in_text', 'zoom_in_tit…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "from ipyleaflet import Map, Marker, Circle\n", + "\n", + "\n", + "MELBOURNE = geolocator.geocode(\"Melbourne\")\n", + "\n", + "\n", + "leaflet_map = Map(center=(MELBOURNE.latitude, MELBOURNE.longitude), zoom=2)\n", + "marker = Marker(location=ORIGIN, draggable=False)\n", + "leaflet_map.add_layer(marker);\n", + "\n", + "for _, coordinates, _, dist, _, _, _, _ in enhanced_data_points:\n", + " marker = Marker(location=coordinates, draggable=False)\n", + " leaflet_map.add_layer(marker);\n", + " \n", + " circle = Circle(location=coordinates, radius=int(dist.meters), fill=False, weight=1)\n", + " leaflet_map.add_layer(circle)\n", + "\n", + "\n", + "display(leaflet_map)\n" + ] + }, + { + "cell_type": "markdown", + "id": "b62287c9", + "metadata": {}, + "source": [ + "However we can fiddle with that so that things end up in roughly the right spot. (Highly scientific. I know...)" + ] + }, + { + "cell_type": "code", + "execution_count": 122, + "id": "1f81d6a2", + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "8f4b34318c5b4d9299ad23184f29b1c0", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Map(center=[-37.8142176, 144.9631608], controls=(ZoomControl(options=['position', 'zoom_in_text', 'zoom_in_tit…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "from ipyleaflet import Map, Marker, Circle\n", + "\n", + "\n", + "FUDGE_FACTOR = 0.821\n", + "\n", + "\n", + "leaflet_map = Map(center=(MELBOURNE.latitude, MELBOURNE.longitude), zoom=2)\n", + "marker = Marker(location=ORIGIN, draggable=False)\n", + "leaflet_map.add_layer(marker);\n", + "\n", + "for _, coordinates, _, dist, _, _, _, _ in enhanced_data_points:\n", + " marker = Marker(location=coordinates, draggable=False)\n", + " leaflet_map.add_layer(marker);\n", + " \n", + " circle = Circle(location=coordinates, radius=int(dist.meters * FUDGE_FACTOR), fill=False, weight=1)\n", + " leaflet_map.add_layer(circle)\n", + "\n", + "\n", + "display(leaflet_map)" + ] + }, + { + "cell_type": "markdown", + "id": "79dd3ef0", + "metadata": {}, + "source": [ + "Also I don't want to mess with the math required to actually propagate back the circles projected on a sphere until they actually intersec.\n", + "\n", + "So let's just assume we have a rough idea when the event occured and we can use that timestamp +-5min to draw some circles and use our eyes to check for intersections." + ] + }, + { + "cell_type": "code", + "execution_count": 135, + "id": "0cb443cc", + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "19841666b00e498e9cf91b0588a12f7c", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Map(center=[-20.5, -175.4], controls=(ZoomControl(options=['position', 'zoom_in_text', 'zoom_in_title', 'zoom_…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "from ipyleaflet import Map, Marker, Circle\n", + "\n", + "\n", + "leaflet_map = Map(center=ORIGIN, zoom=2)\n", + "marker = Marker(location=ORIGIN, draggable=False)\n", + "leaflet_map.add_layer(marker);\n", + "\n", + "\n", + "T_MINUS = datetime(2022, 1, 16, 4, 22, 0, tzinfo=timezone(\"UTC\"))\n", + "T_PLUS = datetime(2022, 1, 16, 4, 32, 0, tzinfo=timezone(\"UTC\"))\n", + "\n", + "\n", + "for location, coordinates, ts, speed in relative_speed_data:\n", + " marker = Marker(location=coordinates, draggable=False)\n", + " leaflet_map.add_layer(marker);\n", + " \n", + " for start, color in [(T_MINUS, \"green\"), (T_ZERO, \"red\"), (T_PLUS, \"blue\")]:\n", + " dist = speed * (ts - start).total_seconds()\n", + "\n", + " circle = Circle(location=coordinates, radius=int(dist * FUDGE_FACTOR), fill=False, weight=1, color=color)\n", + " leaflet_map.add_layer(circle)\n", + "\n", + "\n", + "display(leaflet_map)" + ] + }, + { + "cell_type": "markdown", + "id": "457f9ad9", + "metadata": {}, + "source": [ + "While that looks close, just zoom out until you can see the entire map.\n", + "\n", + "It's a mess.\n", + "Not enough data points, the locations are not spread out wide enough, the timestamps are off, and there is some weird map projection stuff going on...." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5ed3a84b", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.1" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}