Calculate personal count for a bicycle counter
The bicycle counter at the "cycle highway" between Stuttgart-Rohr and Böblingen reached over 1 million cyclists counted. The counting started at 2019-03-30 and I am interested how many times I was counted when cycling there.
First we need the polygon covering the track near the counter (in GeoJSON format):
{"type": "Feature", "properties": {}, "geometry": { "type": "Polygon", "coordinates": [[ [9.081402, 48.711934], [9.082003, 48.711683], [9.081429, 48.711134], [9.080785, 48.711339], [9.081402, 48.711934] ]]}}
The polygon will be loaded into SQLite and used to calculate if my tours intersected with it. To simplify the problem every tour only counts once, even if I crossed the counter more than one time. The following code loads all "fit" files in the folder "tours" (which contains my tours from 2019-05-30 till today); chunks them in lists of 100 dataframes and intersects them with the polygon using Spatialite. The SQLite database is only used to calculate the intersection, so it is only in memory. The libraries required are sqlite-utils, shapely and fitdecode.
import json import warnings from itertools import zip_longest from pathlib import Path from typing import Iterator import fitdecode from shapely.geometry import shape from sqlite_utils.db import Database def get_data(fn: Path) -> Iterator[dict]: with fitdecode.FitReader( fn, processor=fitdecode.StandardUnitsDataProcessor() ) as fit: for frame in fit: if frame.frame_type == fitdecode.FIT_FRAME_DATA: data = {i.name: i.value for i in frame.fields} # only data frames with lat/lon if "position_lat" in data: yield data def process_file(db: Database, fn: Path) -> Iterator[int]: # process 100 dataframes at once to speed up for chunk in zip_longest(*[get_data(fn)] * 100): # get lat/lon and filter empty values coords = [ [c["position_long"], c["position_lat"]] for c in chunk if c and c["position_long"] ] # only use coords with a few values if len(coords) > 2: r = db.query( "select Intersects(GeomFromText(?), Envelope(calc.geometry)) as x from calc", [shape({"type": "LineString", "coordinates": coords}).wkt], ) # [{'x': 1}] -> 1; [{'x': 0}] -> 0 yield list(r)[0]["x"] db = Database(memory=True) db.init_spatialite("/usr/lib/mod_spatialite.so") # path to library file # table to calculate the intersections table = db["calc"] table.insert( { "id": 1, "geojson": { "type": "Polygon", "coordinates": [ [ [9.081402, 48.711934], [9.082003, 48.711683], [9.081429, 48.711134], [9.080785, 48.711339], [9.081402, 48.711934], ] ], }, }, pk="id", ) # add a geometry column and fill with geojson polygon table.add_geometry_column(f"geometry", "POLYGON") for row in table.rows: db.execute( f"update calc set geometry = GeomFromText(?, 4326) where id=?", [shape(json.loads(row["geojson"])).wkt, row["id"]], ) # disable warnings from fitdecode warnings.simplefilter("ignore") counter = 0 files = sorted(Path("tours").glob("*.fit")) for index, fn in enumerate(files, start=1): result = sum(process_file(db, fn)) # only count once, even if sum() > 1 counter += 1 if result else 0 print( f"{index:>4}/{len(files):>4} -- {fn.name:>50} -- {result} -- {counter:>3}" ) print(f"crossed the counter {counter} times")
The last lines of the output:
1579/1581 -- 2023-07-29-081454-ELEMNT BOLT 105D-708-0.fit -- 0 -- 366 1580/1581 -- 2023-07-29-151400-ELEMNT BOLT 105D-709-0.fit -- 1 -- 367 1581/1581 -- 2023-07-30-083422-ELEMNT BOLT 105D-710-0.fit -- 1 -- 368 crossed the counter 368 times
Of the over 1 million cyclists counted, I contributed 368 counts.