diff --git a/resspublica/baselLuftQualitat.py b/resspublica/baselLuftQualitat.py new file mode 100644 index 0000000..59af2f5 --- /dev/null +++ b/resspublica/baselLuftQualitat.py @@ -0,0 +1,410 @@ +from datetime import date +import pandas as pd +import geopandas as gpd +from shapely import wkb +import matplotlib.pyplot as plt + +import logging +logger = logging.getLogger("resspublica") + +from .translations import * +from .utils import * + +station_coordinates = { + "12500": (47.558, 7.588), # https://luftqualitaet.ch/messnetz/station/blSIB + "12450": (47.500, 7.620), # https://luftqualitaet.ch/messnetz/station/blSIB + "12510": (47.450, 7.780), # https://luftqualitaet.ch/messnetz/station/blMUT +} + +def generateBaselLuftqualitat(ASSETS): + logger.info("Generating Luftqualitat map in Basel feed...") + logger.info("Preparing data...") + + urls = [ + "https://data.bs.ch/api/v2/catalog/datasets/100048/exports/parquet", + "https://data.bs.ch/api/v2/catalog/datasets/100050/exports/parquet", + "https://data.bs.ch/api/v2/catalog/datasets/100093/exports/parquet", + "https://data.bs.ch/api/v2/catalog/datasets/100049/exports/parquet", + "https://data.bs.ch/api/v2/catalog/datasets/100178/exports/parquet", + "https://data.bs.ch/api/v2/catalog/datasets/100158/exports/parquet", + "https://data.bl.ch/api/v2/catalog/datasets/12500/exports/parquet", + "https://data.bl.ch/api/v2/catalog/datasets/12450/exports/parquet", + "https://data.bl.ch/api/v2/catalog/datasets/12510/exports/parquet" + ] + + # --------------------------------------------------------- + # Load data + # --------------------------------------------------------- + dataframes = [] + + for url in urls: + logger.debug(f"Querying {url}...") + df = pd.read_parquet(url) + + station_id = url.split("/")[-3] + df["station_id"] = station_id + + # ----------------------------------------------------- + # Normalize datetime to Swiss time + # ----------------------------------------------------- + possible_dates = [ + "datum_zeit", + "timestamp", + "anfangszeit", + "messbeginn" + ] + + date_column = next( + ( + col + for col in possible_dates + if col in df.columns + ), + None + ) + + if date_column is None: + logger.warning( + f"No date column found in {url}, skipping" + ) + continue + + df["date_time"] = pd.to_datetime( + df[date_column], + errors="coerce", + utc=True + ) + + df["date_time"] = ( + df["date_time"] + .dt.tz_convert("Europe/Zurich") + .dt.tz_localize(None) + ) + + + # ----------------------------------------------------- + # Add coordinates for BL datasets + # ----------------------------------------------------- + if station_id in station_coordinates: + + lat, lon = station_coordinates[station_id] + + df["latitude"] = lat + df["longitude"] = lon + + + # ----------------------------------------------------- + # Convert long format datasets + # ----------------------------------------------------- + if ( + "parameter" in df.columns + and "messwert" in df.columns + ): + + df = df.pivot_table( + index=[ + "date_time", + "geo_point_2d" + ], + columns="parameter", + values="messwert", + aggfunc="mean" + ).reset_index() + + + # ----------------------------------------------------- + # Pollutant normalization + # ----------------------------------------------------- + pollutant_mapping = { + + "pm10": [ + "pm10", + "pm10_stundenmittelwerte_ug_m3" + ], + + "pm2_5": [ + "pm2_5", + "pm2.5", + "pm25", + "pm2_5_stundenmittelwerte_ug_m3", + "g107_pm25", + "g125_pm25", + "g131_pm25", + "a2hard_pm25", + "feldbergstr2_pm25", + "stjohann2_pm25" + ], + + "no2": [ + "no2", + "no2_stundenmittelwerte_ug_m3", + "g107_no2", + "g125_no2", + "g131_no2", + "a2hard_no2", + "feldbergstr2_no2", + "stjohann2_no2" + ], + + "o3": [ + "o3", + "o3_stundenmittelwerte_ug_m3", + "g107_o3", + "g107_03", + "g125_o3", + "g131_o3", + "a2hard_o3", + "feldbergstr2_o3", + "stjohann2_o3" + ] + } + + + # ----------------------------------------------------- + # Melt everything into: + # date_time | station_id | pollutant | value | geometry + # ----------------------------------------------------- + parts = [] + + for pollutant, candidates in pollutant_mapping.items(): + + for column in candidates: + + if column not in df.columns: + continue + + keep = [ + "date_time", + "station_id", + column + ] + + for extra in [ + "geo_point_2d", + "latitude", + "longitude" + ]: + if extra in df.columns: + keep.append(extra) + + + tmp = df[keep].copy() + + tmp = tmp.rename( + columns={ + column: "value" + } + ) + + tmp["pollutant"] = pollutant + + parts.append(tmp) + + + if not parts: + logger.warning( + f"No pollutants found in {url}" + ) + continue + + + df = pd.concat( + parts, + ignore_index=True + ) + + dataframes.append(df) + + + # --------------------------------------------------------- + # Combine datasets + # --------------------------------------------------------- + dataframe = pd.concat( + dataframes, + ignore_index=True, + sort=False + ) + + + # --------------------------------------------------------- + # Geometry + # --------------------------------------------------------- + def safe_load(x): + try: + return wkb.loads(x) + except Exception: + return None + + + dataframe["geometry"] = None + + + if "geo_point_2d" in dataframe.columns: + + dataframe["geometry"] = dataframe[ + "geo_point_2d" + ].apply( + safe_load + ) + + + # Fill missing geometry from coordinates + missing_geometry = dataframe["geometry"].isna() + + dataframe.loc[ + missing_geometry, + "geometry" + ] = gpd.points_from_xy( + dataframe.loc[missing_geometry, "longitude"], + dataframe.loc[missing_geometry, "latitude"] + ) + + + geo = gpd.GeoDataFrame( + dataframe, + geometry="geometry", + crs="EPSG:4326" + ) + + + # --------------------------------------------------------- + # Time filter + # --------------------------------------------------------- + start = pd.Timestamp( + "2026-07-29 00:00:00" + ) + + end = pd.Timestamp( + "2026-07-29 23:59:59" + ) + + + geo = geo[ + (geo["date_time"] >= start) + & + (geo["date_time"] <= end) + ] + + + # --------------------------------------------------------- + # Average per station + # --------------------------------------------------------- + averaged = ( + geo + .groupby( + [ + "station_id", + "pollutant", + "geometry" + ], + as_index=False + ) + ["value"] + .mean() + ) + + + averaged = gpd.GeoDataFrame( + averaged, + geometry="geometry", + crs="EPSG:4326" + ) + # --------------------------------------------------------- + # Boundaries + # --------------------------------------------------------- + gdb = ( + ASSETS / + "swissBOUNDARIES3D_1_5_LV95_LN02.gdb" + ) + + cantons = gpd.read_file( + gdb, + layer="TLM_KANTONSGEBIET" + ) + + basel_stadt = cantons[ + cantons["KANTONSNUMMER"] == 12 + ].to_crs(averaged.crs) + + basel_land = cantons[ + cantons["KANTONSNUMMER"] == 13 + ].to_crs(averaged.crs) + + # --------------------------------------------------------- + # Plot + # --------------------------------------------------------- + pollutants = [ + "pm10", + "pm2_5", + "no2", + "o3" + ] + + + fig, axes = plt.subplots( + 1, + len(pollutants), + figsize=(20, 5) + ) + + + for ax, pollutant in zip( + axes, + pollutants + ): + + subset = averaged[ + averaged["pollutant"] == pollutant + ] + + + if subset.empty: + ax.set_visible(False) + continue + + + subset.plot( + ax=ax, + column="value", + cmap="hot", + legend=True, + markersize=80 + ) + + + basel_land.plot( + ax=ax, + facecolor="none", + edgecolor="black", + linewidth=2 + ) + + + basel_stadt.plot( + ax=ax, + facecolor="none", + edgecolor="black", + linewidth=1 + ) + + + for _, row in subset.iterrows(): + ax.annotate( + row["station_id"], + ( + row.geometry.x, + row.geometry.y + ), + fontsize=8 + ) + + + ax.set_title( + pollutant + ) + + ax.axis("off") + + + plt.tight_layout() + plt.show() diff --git a/resspublica/main.py b/resspublica/main.py index af03317..db5898f 100644 --- a/resspublica/main.py +++ b/resspublica/main.py @@ -10,12 +10,14 @@ from .federalInitiativesFeeds import * from .translations import * from .bernAsianHornet import * from .bernReligionMap import * +from .baselLuftQualitat import * def main(): parser = argparse.ArgumentParser() parser.add_argument("--verbose", action="store_true", help="Enable debug logging") parser.add_argument("--gen_federalInitiatives", action="store_true", help="Generate federal popular initiatives feed") parser.add_argument("--gen_bernAsianHornets", action="store_true", help="Generate Asian hornets sightings in Bern feed (only Mondays)") + parser.add_argument("--gen_baselLuftqualitat", action="store_true", help="Generate Basel air quality feed") parser.add_argument("--force_gen_bernAsianHornets", action="store_true", help="Generate Asian hornets sightings in Bern feed even when not Monday") parser.add_argument("--gen_bernReligionMap", action="store_true", help="Generate map of religions in Bern (only start of trimesters)") parser.add_argument("--force_gen_bernReligionMap", action="store_true", help="Generate map of religions in Bern even when not start of trimester") @@ -40,4 +42,6 @@ def main(): generateBernAsianHornetFeed(ASSETS, CACHE) if ( args.gen_bernReligionMap and date.today().day == 1 and date.today().month in [1, 5, 9]) or args.force_gen_bernReligionMap: generateBernReligionMap(ASSETS, CACHE) + if args.gen_baselLuftqualitat: + generateBaselLuftqualitat(ASSETS) logging.info("Done")